Magnetic_Field.cc

Go to the documentation of this file.
00001 //! \file
00002 //!
00003 //! Source file for Magnetic_Field class.
00004 //!
00005 //! This defines the Magnetic_Field class and the member routines which
00006 //! construct the OLYMPUS toroidal magnetic field.
00007 //!
00008 //! \author D.K. Hasell
00009 //! \version 1.0
00010 //! \date 2010-10-14
00011 //!
00012 //! \ingroup detector
00013 
00014 // Include the Magnetic_Field and other user header files.
00015 
00016 #include "Magnetic_Field.h"
00017 
00018 #include "MF_Messenger.h"
00019 
00020 // Include the GEANT4 header files referenced here.
00021 
00022 #include "globals.hh"
00023 #include "G4ios.hh"
00024 
00025 // Include the C++ header files referenced here.
00026 
00027 #include <fstream>
00028 
00029 // Constructor for Magnetic_Field class.
00030 
00031 Magnetic_Field::Magnetic_Field( const char * filename, G4double scale )
00032    : fScale( scale ) {
00033 
00034    // Pass pointer to this class to Messenger.
00035 
00036    MF_Messenger::Instance()->setMFptr( this );
00037 
00038    G4cout << "     reading field grid\n" << flush;
00039 
00040    // Open file with grid data.
00041 
00042    ifstream grid( filename );
00043 
00044    // Read the number of steps, starting value, and step size for grid.
00045 
00046    grid >> fNdata >> fNx >> fNy >> fNz
00047         >> fXmin >> fYmin >> fZmin
00048         >> fDx >> fDy >> fDz;
00049 
00050    // Convert to proper units.
00051 
00052    fXmin *= millimeter;
00053    fYmin *= millimeter;
00054    fZmin *= millimeter;
00055 
00056    fDx *= millimeter;
00057    fDy *= millimeter;
00058    fDz *= millimeter;
00059 
00060    // Resize the vectors which store the grid data.
00061 
00062    fBx.resize( fNx );
00063    fBy.resize( fNx );
00064    fBz.resize( fNx );
00065 
00066    for( G4int ix = 0; ix < fNx; ++ix ){
00067 
00068       fBx[ix].resize( fNy );
00069       fBy[ix].resize( fNy );
00070       fBz[ix].resize( fNy );
00071 
00072       for( G4int iy = 0; iy < fNy; ++iy ) {
00073 
00074          fBx[ix][iy].resize( fNz );
00075          fBy[ix][iy].resize( fNz );
00076          fBz[ix][iy].resize( fNz );
00077 
00078       }
00079 
00080    }
00081 
00082    // Read field data and fill vectors with field components in proper units.
00083 
00084    double bx, by, bz;
00085 
00086    for( G4int iz = 0; iz < fNz; ++iz ) {
00087       for( G4int iy = 0; iy < fNy; ++iy ) {
00088          for( G4int ix = 0; ix < fNx; ++ix ) {
00089 
00090             grid >> bx >> by >> bz;
00091 
00092             fBx[ix][iy][iz] = bx * gauss;
00093             fBy[ix][iy][iz] = by * gauss;
00094             fBz[ix][iy][iz] = bz * gauss;
00095 
00096          }
00097       }
00098    }
00099 
00100    // Close data file.
00101 
00102    grid.close();
00103 
00104 }
00105 
00106 // Member function to return the magnetic field value at a given point.
00107 
00108 void Magnetic_Field::GetFieldValue( const G4double Point[4],
00109                                     G4double * Bfield ) const {
00110 
00111    G4double x = Point[0];
00112    G4double y = Point[1];
00113    G4double z = Point[2];
00114 
00115    // Find grid cubic containing the point.
00116 
00117    G4int ix0 = ( x - fXmin ) / fDx;
00118    G4int iy0 = ( y - fYmin ) / fDy;
00119    G4int iz0 = ( z - fZmin ) / fDz;
00120 
00121    G4int ix1 = ix0 + 1;
00122    G4int iy1 = iy0 + 1;
00123    G4int iz1 = iz0 + 1;
00124 
00125    // Declare fractions of interval.
00126 
00127    G4double fx0, fx1, fy0, fy1, fz0, fz1;
00128 
00129    // Test that point is within the range of the grid.
00130 
00131    if( ix1 > 0 && ix1 < fNx &&
00132        iy1 > 0 && iy1 < fNy &&
00133        iz1 > 0 && iz1 < fNz ) {
00134 
00135       // Calculate the fractions of the interval for the given point.
00136 
00137       fx0 = ( x - fXmin - ix0 * fDx ) / fDx;
00138       fy0 = ( y - fYmin - iy0 * fDy ) / fDy;
00139       fz0 = ( z - fZmin - iz0 * fDz ) / fDz;
00140 
00141       fx1 = ( fXmin + ix1 * fDx - x ) / fDx;
00142       fy1 = ( fYmin + iy1 * fDy - y ) / fDy;
00143       fz1 = ( fZmin + iz1 * fDz - z ) / fDz;
00144 
00145       // Interpolate the magnetic field from the values at the 8 corners.
00146 
00147       Bfield[0] = fz1*(fy1*( fx1*fBx[ix0][iy0][iz0]+fx0*fBx[ix1][iy0][iz0] )
00148                        + fy0*( fx1*fBx[ix0][iy1][iz0]+fx0*fBx[ix1][iy1][iz0] ))
00149          + fz0*(fy1*( fx1*fBx[ix0][iy0][iz1]+fx0*fBx[ix1][iy0][iz1] )
00150                 + fy0*( fx1*fBx[ix0][iy1][iz1]+fx0*fBx[ix1][iy1][iz1] ));
00151       
00152       Bfield[1] = fz1*(fy1*( fx1*fBy[ix0][iy0][iz0]+fx0*fBy[ix1][iy0][iz0] )
00153                        + fy0*( fx1*fBy[ix0][iy1][iz0]+fx0*fBy[ix1][iy1][iz0] ))
00154          + fz0*(fy1*( fx1*fBy[ix0][iy0][iz1]+fx0*fBy[ix1][iy0][iz1] )
00155                 + fy0*( fx1*fBy[ix0][iy1][iz1]+fx0*fBy[ix1][iy1][iz1] ));
00156 
00157       Bfield[2] = fz1*(fy1*( fx1*fBz[ix0][iy0][iz0]+fx0*fBz[ix1][iy0][iz0] )
00158                        + fy0*( fx1*fBz[ix0][iy1][iz0]+fx0*fBz[ix1][iy1][iz0] ))
00159          + fz0*(fy1*( fx1*fBz[ix0][iy0][iz1]+fx0*fBz[ix1][iy0][iz1] )
00160                 + fy0*( fx1*fBz[ix0][iy1][iz1]+fx0*fBz[ix1][iy1][iz1] ));
00161 
00162    }
00163 
00164    else {
00165       Bfield[0] = 0.0;
00166       Bfield[1] = 0.0;
00167       Bfield[2] = 0.0;
00168    }
00169 
00170    // Scale field value.
00171 
00172    Bfield[0] *= fScale;
00173    Bfield[1] *= fScale;
00174    Bfield[2] *= fScale;
00175 
00176    // This is a pure magnetic field so set electric field to zero just in case.
00177 
00178    Bfield[3] = 0.0;
00179    Bfield[4] = 0.0;
00180    Bfield[5] = 0.0;
00181 
00182    return;
00183 
00184 }
00185 
00186 // Set the magnetic field scaling factor.
00187 
00188 G4double Magnetic_Field::setScale( G4double scale ) {
00189    return fScale = scale;
00190 }
00191 
00192 // Get the magnetic field scaling factor.
00193 
00194 G4double Magnetic_Field::getScale() {
00195    return fScale;
00196 }