Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/hadrontherapy/src/HadrontherapyMagneticField3D.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 ]

Diff markup

Differences between /examples/advanced/hadrontherapy/src/HadrontherapyMagneticField3D.cc (Version 11.3.0) and /examples/advanced/hadrontherapy/src/HadrontherapyMagneticField3D.cc (Version 9.1.p2)


  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 "HadrontherapyMagneticField3D.hh"        
 30 #include "G4SystemOfUnits.hh"                     
 31 #include "G4AutoLock.hh"                          
 32                                                   
 33 namespace{  G4Mutex MyHadrontherapyLock=G4MUTE    
 34                                                   
 35 using namespace std;                              
 36                                                   
 37 HadrontherapyMagneticField3D::HadrontherapyMag    
 38   :fXoffset(xOffset),invertX(false),invertY(fa    
 39 {                                                 
 40    //The format file is: X Y Z Ex Ey Ez           
 41                                                   
 42   double lenUnit= meter;                          
 43   double fieldUnit= tesla;                        
 44   G4cout << "\n-------------------------------    
 45    << "\n      Magnetic field"                    
 46    << "\n-------------------------------------    
 47                                                   
 48                                                   
 49   G4cout << "\n ---> " "Reading the field grid    
 50   G4AutoLock lock(&MyHadrontherapyLock);          
 51                                                   
 52   ifstream file( filename ); // Open the file     
 53                                                   
 54   // Ignore first blank line                      
 55   char buffer[256];                               
 56   file.getline(buffer,256);                       
 57                                                   
 58   // Read table dimensions                        
 59   file >> nx >> ny >> nz; // Note dodgy order     
 60                                                   
 61   G4cout << "  [ Number of values x,y,z: "        
 62    << nx << " " << ny << " " << nz << " ] "       
 63    << G4endl;                                     
 64                                                   
 65   // Set up storage space for table               
 66   xField.resize( nx );                            
 67   yField.resize( nx );                            
 68   zField.resize( nx );                            
 69   int ix, iy, iz;                                 
 70   for (ix=0; ix<nx; ix++) {                       
 71     xField[ix].resize(ny);                        
 72     yField[ix].resize(ny);                        
 73     zField[ix].resize(ny);                        
 74     for (iy=0; iy<ny; iy++) {                     
 75       xField[ix][iy].resize(nz);                  
 76       yField[ix][iy].resize(nz);                  
 77       zField[ix][iy].resize(nz);                  
 78     }                                             
 79   }                                               
 80                                                   
 81   // Read in the data                             
 82   G4double xval=0.;                               
 83   G4double yval=0.;                               
 84   G4double zval=0.;                               
 85   G4double bx=0.;                                 
 86   G4double by=0.;                                 
 87   G4double bz=0.;                                 
 88   for (ix=0; ix<nx; ix++) {                       
 89     for (iy=0; iy<ny; iy++) {                     
 90       for (iz=0; iz<nz; iz++) {                   
 91         file >> xval >> yval >> zval >> bx >>     
 92         if ( ix==0 && iy==0 && iz==0 ) {          
 93           minx = xval * lenUnit;                  
 94           miny = yval * lenUnit;                  
 95           minz = zval * lenUnit;                  
 96         }                                         
 97         xField[ix][iy][iz] = bx * fieldUnit;      
 98         yField[ix][iy][iz] = by * fieldUnit;      
 99         zField[ix][iy][iz] = bz * fieldUnit;      
100       }                                           
101     }                                             
102   }                                               
103   file.close();                                   
104                                                   
105   lock.unlock();                                  
106                                                   
107   maxx = xval * lenUnit;                          
108   maxy = yval * lenUnit;                          
109   maxz = zval * lenUnit;                          
110                                                   
111   G4cout << "\n ---> ... done reading " << G4e    
112                                                   
113   // G4cout << " Read values of field from fil    
114   G4cout << " ---> assumed the order:  x, y, z    
115    << "\n ---> Min values x,y,z: "                
116    << minx/cm << " " << miny/cm << " " << minz    
117    << "\n ---> Max values x,y,z: "                
118    << maxx/cm << " " << maxy/cm << " " << maxz    
119    << "\n ---> The field will be offset by " <    
120                                                   
121   // Should really check that the limits are n    
122   if (maxx < minx) {swap(maxx,minx); invertX =    
123   if (maxy < miny) {swap(maxy,miny); invertY =    
124   if (maxz < minz) {swap(maxz,minz); invertZ =    
125   G4cout << "\nAfter reordering if neccesary"     
126    << "\n ---> Min values x,y,z: "                
127    << minx/cm << " " << miny/cm << " " << minz    
128    << " \n ---> Max values x,y,z: "               
129    << maxx/cm << " " << maxy/cm << " " << maxz    
130                                                   
131   dx = maxx - minx;                               
132   dy = maxy - miny;                               
133   dz = maxz - minz;                               
134   G4cout << "\n ---> Dif values x,y,z (range):    
135    << dx/cm << " " << dy/cm << " " << dz/cm <<    
136    << "\n-------------------------------------    
137 }                                                 
138                                                   
139 void HadrontherapyMagneticField3D::GetFieldVal    
140               double *Bfield ) const              
141 {                                                 
142     double x = point[0]+ fXoffset;                
143     double y = point[1];                          
144     double z = point[2];                          
145                                                   
146     // Position of given point within region,     
147     // [0,1]                                      
148     double xfraction = (x - minx) / dx;           
149     double yfraction = (y - miny) / dy;           
150     double zfraction = (z - minz) / dz;           
151                                                   
152     if (invertX) { xfraction = 1 - xfraction;}    
153     if (invertY) { yfraction = 1 - yfraction;}    
154     if (invertZ) { zfraction = 1 - zfraction;}    
155                                                   
156     // Need addresses of these to pass to modf    
157     // modf uses its second argument as an OUT    
158     double xdindex, ydindex, zdindex;             
159                                                   
160     // Position of the point within the cuboid    
161     // nearest surrounding tabulated points       
162     double xlocal = ( std::modf(xfraction*(nx-    
163     double ylocal = ( std::modf(yfraction*(ny-    
164     double zlocal = ( std::modf(zfraction*(nz-    
165                                                   
166     // The indices of the nearest tabulated po    
167     // are all less than those of the given po    
168     int xindex = static_cast<int>(std::floor(x    
169     int yindex = static_cast<int>(std::floor(y    
170     int zindex = static_cast<int>(std::floor(z    
171                                                   
172       // Check that the point is within the de    
173     if ((xindex < 0) || (xindex >= nx - 1) ||     
174         (yindex < 0) || (yindex >= ny - 1) ||     
175         (zindex < 0) || (zindex >= nz - 1))       
176     {                                             
177         Bfield[0] = 0.0;                          
178         Bfield[1] = 0.0;                          
179         Bfield[2] = 0.0;                          
180     }                                             
181     else                                          
182     {                                             
183                                                   
184 #ifdef DEBUG_INTERPOLATING_FIELD                  
185         G4cout << "Local x,y,z: " << xlocal <<    
186         G4cout << "Index x,y,z: " << xindex <<    
187         double valx0z0, mulx0z0, valx1z0, mulx    
188         double valx0z1, mulx0z1, valx1z1, mulx    
189         valx0z0= table[xindex  ][0][zindex];      
190         valx1z0= table[xindex+1][0][zindex];      
191         valx0z1= table[xindex  ][0][zindex+1];    
192         valx1z1= table[xindex+1][0][zindex+1];    
193 #endif                                            
194                                                   
195         // Full 3-dimensional version             
196         Bfield[0] =                               
197           xField[xindex  ][yindex  ][zindex  ]    
198           xField[xindex  ][yindex  ][zindex+1]    
199           xField[xindex  ][yindex+1][zindex  ]    
200           xField[xindex  ][yindex+1][zindex+1]    
201           xField[xindex+1][yindex  ][zindex  ]    
202           xField[xindex+1][yindex  ][zindex+1]    
203           xField[xindex+1][yindex+1][zindex  ]    
204           xField[xindex+1][yindex+1][zindex+1]    
205                                                   
206         Bfield[1] =                               
207           yField[xindex  ][yindex  ][zindex  ]    
208           yField[xindex  ][yindex  ][zindex+1]    
209           yField[xindex  ][yindex+1][zindex  ]    
210           yField[xindex  ][yindex+1][zindex+1]    
211           yField[xindex+1][yindex  ][zindex  ]    
212           yField[xindex+1][yindex  ][zindex+1]    
213           yField[xindex+1][yindex+1][zindex  ]    
214           yField[xindex+1][yindex+1][zindex+1]    
215                                                   
216         Bfield[2] =                               
217           zField[xindex  ][yindex  ][zindex  ]    
218           zField[xindex  ][yindex  ][zindex+1]    
219           zField[xindex  ][yindex+1][zindex  ]    
220           zField[xindex  ][yindex+1][zindex+1]    
221           zField[xindex+1][yindex  ][zindex  ]    
222           zField[xindex+1][yindex  ][zindex+1]    
223           zField[xindex+1][yindex+1][zindex  ]    
224           zField[xindex+1][yindex+1][zindex+1]    
225     }                                             
226                                                   
227 //In order to obtain the output file with the     
228 /*  std::ofstream MagneticField("MagneticField    
229      MagneticField<<   Bfield[0] << '\t' << "     
230       <<   Bfield[1] << '\t' << "    "            
231       <<   Bfield[2] << '\t' << "   "             
232       <<   point[0] << '\t' << "   "              
233       <<   point[1] << '\t' << "    "             
234       <<   point[2] << '\t' << "   "              
235       << std::endl;*/                             
236                                                   
237 }                                                 
238