#include <HadrontherapyMagneticField3D.hh>

Inheritance diagram for HadrontherapyMagneticField3D:

Collaboration diagram for HadrontherapyMagneticField3D:

Public Member Functions
	HadrontherapyMagneticField3D (const char *filename, double xOffset)

void	GetFieldValue (const double Point[4], double *Bfield) const

Public Member Functions inherited from G4MagneticField
	G4MagneticField ()

virtual	~G4MagneticField ()

	G4MagneticField (const G4MagneticField &r)

G4MagneticField &	operator= (const G4MagneticField &p)

G4bool	DoesFieldChangeEnergy () const

Public Member Functions inherited from G4ElectroMagneticField
	G4ElectroMagneticField ()

virtual	~G4ElectroMagneticField ()

	G4ElectroMagneticField (const G4ElectroMagneticField &r)

G4ElectroMagneticField &	operator= (const G4ElectroMagneticField &p)

Public Member Functions inherited from G4Field
	G4Field (G4bool gravityOn=false)

	G4Field (const G4Field &)

virtual	~G4Field ()

G4Field &	operator= (const G4Field &p)

G4bool	IsGravityActive () const

void	SetGravityActive (G4bool OnOffFlag)

virtual G4Field *	Clone () const

Private Attributes
vector< vector< vector< double > > >	xField

vector< vector< vector< double > > >	yField

vector< vector< vector< double > > >	zField

int	nx

int	ny

int	nz

double	minx

double	maxx

double	miny

double	maxy

double	minz

double	maxz

double	dx

double	dy

double	dz

double	fXoffset

bool	invertX

bool	invertY

bool	invertZ

Detailed Description

Definition at line 40 of file HadrontherapyMagneticField3D.hh.

Constructor & Destructor Documentation

◆ HadrontherapyMagneticField3D()

HadrontherapyMagneticField3D::HadrontherapyMagneticField3D	(	const char *	filename,
		double	xOffset
	)

Definition at line 35 of file HadrontherapyMagneticField3D.cc.

   :fXoffset(xOffset),invertX(false),invertY(false),invertZ(false)
 {    
    //The format file is: X Y Z Ex Ey Ez
 
   double lenUnit= meter; 
   double fieldUnit= tesla;  
   G4cout << "\n-----------------------------------------------------------"
      << "\n      Magnetic field"
      << "\n-----------------------------------------------------------";
 
 
   G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << endl; 
   G4AutoLock lock(&MyHadrontherapyLock);
 
   ifstream file( filename ); // Open the file for reading.
   
   // Ignore first blank line
   char buffer[256];
   file.getline(buffer,256);
   
   // Read table dimensions 
   file >> nx >> ny >> nz; // Note dodgy order
 
   G4cout << "  [ Number of values x,y,z: " 
      << nx << " " << ny << " " << nz << " ] "
      << endl;
 
   // Set up storage space for table
   xField.resize( nx );
   yField.resize( nx );
   zField.resize( nx );
   int ix, iy, iz;
   for (ix=0; ix<nx; ix++) {
     xField[ix].resize(ny);
     yField[ix].resize(ny);
     zField[ix].resize(ny);
     for (iy=0; iy<ny; iy++) {
       xField[ix][iy].resize(nz);
       yField[ix][iy].resize(nz);
       zField[ix][iy].resize(nz);
     }
   }
 
   // Read in the data
   G4double xval=0.;
   G4double yval=0.;
   G4double zval=0.;
   G4double bx=0.;
   G4double by=0.;
   G4double bz=0.;
   for (ix=0; ix<nx; ix++) {
     for (iy=0; iy<ny; iy++) {
       for (iz=0; iz<nz; iz++) {
         file >> xval >> yval >> zval >> bx >> by >> bz ;
         if ( ix==0 && iy==0 && iz==0 ) {
           minx = xval * lenUnit;
           miny = yval * lenUnit;
           minz = zval * lenUnit;
         }
         xField[ix][iy][iz] = bx * fieldUnit; 
         yField[ix][iy][iz] = by * fieldUnit;
         zField[ix][iy][iz] = bz * fieldUnit;
       }
     }
   }
   file.close();
 
   lock.unlock();
 
   maxx = xval * lenUnit;
   maxy = yval * lenUnit;
   maxz = zval * lenUnit;
 
   G4cout << "\n ---> ... done reading " << endl;
 
   // G4cout << " Read values of field from file " << filename << endl; 
   G4cout << " ---> assumed the order:  x, y, z, Bx, By, Bz "
      << "\n ---> Min values x,y,z: " 
      << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
      << "\n ---> Max values x,y,z: " 
      << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm "
      << "\n ---> The field will be offset by " << xOffset/cm << " cm " << endl;
 
   // Should really check that the limits are not the wrong way around.
   if (maxx < minx) {swap(maxx,minx); invertX = true;} 
   if (maxy < miny) {swap(maxy,miny); invertY = true;} 
   if (maxz < minz) {swap(maxz,minz); invertZ = true;} 
   G4cout << "\nAfter reordering if neccesary"  
      << "\n ---> Min values x,y,z: " 
      << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
      << " \n ---> Max values x,y,z: " 
      << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm ";
 
   dx = maxx - minx;
   dy = maxy - miny;
   dz = maxz - minz;
   G4cout << "\n ---> Dif values x,y,z (range): " 
      << dx/cm << " " << dy/cm << " " << dz/cm << " cm in z "
      << "\n-----------------------------------------------------------" << endl;
 }

Here is the call graph for this function:

Member Function Documentation

◆ GetFieldValue()

void HadrontherapyMagneticField3D::GetFieldValue	(	const double	Point[4],
		double *	Bfield
	)		const

virtual

Implements G4MagneticField.

Definition at line 137 of file HadrontherapyMagneticField3D.cc.

 {
   double x = point[0]+ fXoffset;
   double y = point[1];
   double z = point[2];
 
   // Check that the point is within the defined region 
   if ( x>=minx && x<maxx &&
        y>=miny && y<maxy && 
        z>=minz && z<maxz ) {
     // Position of given point within region, normalized to the range
     // [0,1]
     double xfraction = (x - minx) / dx;
     double yfraction = (y - miny) / dy; 
     double zfraction = (z - minz) / dz;
 
     if (invertX) { xfraction = 1 - xfraction;}
     if (invertY) { yfraction = 1 - yfraction;}
     if (invertZ) { zfraction = 1 - zfraction;}
 
     // Need addresses of these to pass to modf below.
     // modf uses its second argument as an OUTPUT argument.
     double xdindex, ydindex, zdindex;
     
     // Position of the point within the cuboid defined by the
     // nearest surrounding tabulated points
     double xlocal = ( std::modf(xfraction*(nx-1), &xdindex));
     double ylocal = ( std::modf(yfraction*(ny-1), &ydindex));
     double zlocal = ( std::modf(zfraction*(nz-1), &zdindex));
     
     // The indices of the nearest tabulated point whose coordinates
     // are all less than those of the given point
     int xindex = static_cast<int>(xdindex);
     int yindex = static_cast<int>(ydindex);
     int zindex = static_cast<int>(zdindex);
     
 
 #ifdef DEBUG_INTERPOLATING_FIELD
     G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << endl;
     G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << endl;
     double valx0z0, mulx0z0, valx1z0, mulx1z0;
     double valx0z1, mulx0z1, valx1z1, mulx1z1;
     valx0z0= table[xindex  ][0][zindex];  mulx0z0=  (1-xlocal) * (1-zlocal);
     valx1z0= table[xindex+1][0][zindex];  mulx1z0=   xlocal    * (1-zlocal);
     valx0z1= table[xindex  ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
     valx1z1= table[xindex+1][0][zindex+1]; mulx1z1=  xlocal    * zlocal;
 #endif
         // Full 3-dimensional version
     Bfield[0] =
       xField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
       xField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
       xField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
       xField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
       xField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
       xField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
       xField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
       xField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
 
     Bfield[1] =
       yField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
       yField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
       yField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
       yField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
       yField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
       yField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
       yField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
       yField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
 
     Bfield[2] =
       zField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
       zField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
       zField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
       zField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
       zField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
       zField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
       zField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
       zField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
     
   } else {
     Bfield[0] = 0.0;
     Bfield[1] = 0.0;
     Bfield[2] = 0.0;
   }
 //In order to obtain the output file with the magnetic components read from a particle passing in the magnetic field
 /*  std::ofstream MagneticField("MagneticField.out", std::ios::app);
        MagneticField<<   Bfield[0] << '\t' << "   "
             <<   Bfield[1] << '\t' << "    "
             <<   Bfield[2] << '\t' << "   "
             <<   point[0] << '\t' << "   "
             <<   point[1] << '\t' << "    "
             <<   point[2] << '\t' << "   "
             << G4endl;*/
     
 }