Magnetic_Field.cc

Go to the documentation of this file.

//! \file
//!
//! Source file for Magnetic_Field class.
//!
//! This defines the Magnetic_Field class and the member routines which
//! construct the OLYMPUS toroidal magnetic field.
//!
//! \author D.K. Hasell
//! \version 1.0
//! \date 2010-10-14
//!
//! \ingroup detector

// Include the Magnetic_Field and other user header files.

#include "Magnetic_Field.h"

#include "MF_Messenger.h"

// Include the GEANT4 header files referenced here.

#include "globals.hh"
#include "G4ios.hh"

// Include the C++ header files referenced here.

#include <fstream>

// Constructor for Magnetic_Field class.

Magnetic_Field::Magnetic_Field( const char * filename, G4double scale )
   : fScale( scale ) {

   // Pass pointer to this class to Messenger.

   MF_Messenger::Instance()->setMFptr( this );

   G4cout << "     reading field grid\n" << flush;

   // Open file with grid data.

   ifstream grid( filename );

   // Read the number of steps, starting value, and step size for grid.

   grid >> fNdata >> fNx >> fNy >> fNz
        >> fXmin >> fYmin >> fZmin
        >> fDx >> fDy >> fDz;

   // Convert to proper units.

   fXmin *= millimeter;
   fYmin *= millimeter;
   fZmin *= millimeter;

   fDx *= millimeter;
   fDy *= millimeter;
   fDz *= millimeter;

   // Resize the vectors which store the grid data.

   fBx.resize( fNx );
   fBy.resize( fNx );
   fBz.resize( fNx );

   for( G4int ix = 0; ix < fNx; ++ix ){

      fBx[ix].resize( fNy );
      fBy[ix].resize( fNy );
      fBz[ix].resize( fNy );

      for( G4int iy = 0; iy < fNy; ++iy ) {

         fBx[ix][iy].resize( fNz );
         fBy[ix][iy].resize( fNz );
         fBz[ix][iy].resize( fNz );

      }

   }

   // Read field data and fill vectors with field components in proper units.

   double bx, by, bz;

   for( G4int iz = 0; iz < fNz; ++iz ) {
      for( G4int iy = 0; iy < fNy; ++iy ) {
         for( G4int ix = 0; ix < fNx; ++ix ) {

            grid >> bx >> by >> bz;

            fBx[ix][iy][iz] = bx * gauss;
            fBy[ix][iy][iz] = by * gauss;
            fBz[ix][iy][iz] = bz * gauss;

         }
      }
   }

   // Close data file.

   grid.close();

}

// Member function to return the magnetic field value at a given point.

void Magnetic_Field::GetFieldValue( const G4double Point[4],
                                    G4double * Bfield ) const {

   G4double x = Point[0];
   G4double y = Point[1];
   G4double z = Point[2];

   // Find grid cubic containing the point.

   G4int ix0 = ( x - fXmin ) / fDx;
   G4int iy0 = ( y - fYmin ) / fDy;
   G4int iz0 = ( z - fZmin ) / fDz;

   G4int ix1 = ix0 + 1;
   G4int iy1 = iy0 + 1;
   G4int iz1 = iz0 + 1;

   // Declare fractions of interval.

   G4double fx0, fx1, fy0, fy1, fz0, fz1;

   // Test that point is within the range of the grid.

   if( ix1 > 0 && ix1 < fNx &&
       iy1 > 0 && iy1 < fNy &&
       iz1 > 0 && iz1 < fNz ) {

      // Calculate the fractions of the interval for the given point.

      fx0 = ( x - fXmin - ix0 * fDx ) / fDx;
      fy0 = ( y - fYmin - iy0 * fDy ) / fDy;
      fz0 = ( z - fZmin - iz0 * fDz ) / fDz;

      fx1 = ( fXmin + ix1 * fDx - x ) / fDx;
      fy1 = ( fYmin + iy1 * fDy - y ) / fDy;
      fz1 = ( fZmin + iz1 * fDz - z ) / fDz;

      // Interpolate the magnetic field from the values at the 8 corners.

      Bfield[0] = fz1*(fy1*( fx1*fBx[ix0][iy0][iz0]+fx0*fBx[ix1][iy0][iz0] )
                       + fy0*( fx1*fBx[ix0][iy1][iz0]+fx0*fBx[ix1][iy1][iz0] ))
         + fz0*(fy1*( fx1*fBx[ix0][iy0][iz1]+fx0*fBx[ix1][iy0][iz1] )
                + fy0*( fx1*fBx[ix0][iy1][iz1]+fx0*fBx[ix1][iy1][iz1] ));
      
      Bfield[1] = fz1*(fy1*( fx1*fBy[ix0][iy0][iz0]+fx0*fBy[ix1][iy0][iz0] )
                       + fy0*( fx1*fBy[ix0][iy1][iz0]+fx0*fBy[ix1][iy1][iz0] ))
         + fz0*(fy1*( fx1*fBy[ix0][iy0][iz1]+fx0*fBy[ix1][iy0][iz1] )
                + fy0*( fx1*fBy[ix0][iy1][iz1]+fx0*fBy[ix1][iy1][iz1] ));

      Bfield[2] = fz1*(fy1*( fx1*fBz[ix0][iy0][iz0]+fx0*fBz[ix1][iy0][iz0] )
                       + fy0*( fx1*fBz[ix0][iy1][iz0]+fx0*fBz[ix1][iy1][iz0] ))
         + fz0*(fy1*( fx1*fBz[ix0][iy0][iz1]+fx0*fBz[ix1][iy0][iz1] )
                + fy0*( fx1*fBz[ix0][iy1][iz1]+fx0*fBz[ix1][iy1][iz1] ));

   }

   else {
      Bfield[0] = 0.0;
      Bfield[1] = 0.0;
      Bfield[2] = 0.0;
   }

   // Scale field value.

   Bfield[0] *= fScale;
   Bfield[1] *= fScale;
   Bfield[2] *= fScale;

   // This is a pure magnetic field so set electric field to zero just in case.

   Bfield[3] = 0.0;
   Bfield[4] = 0.0;
   Bfield[5] = 0.0;

   return;

}

// Set the magnetic field scaling factor.

G4double Magnetic_Field::setScale( G4double scale ) {
   return fScale = scale;
}

// Get the magnetic field scaling factor.

G4double Magnetic_Field::getScale() {
   return fScale;
}