00001
00002 #include <iostream>
00003
00004 #include "BDSGlobalConstants.hh"
00005 #include "BDS3DMagField.hh"
00006
00007 BDS3DMagField::BDS3DMagField( const char* filename, double zOffset )
00008 :fZoffset(zOffset),invertX(false),invertY(false),invertZ(false)
00009 {
00010 _lenUnit= CLHEP::cm;
00011 _fieldUnit= CLHEP::tesla;
00012
00013 G4cout << "\n-----------------------------------------------------------"
00014 << "\n Magnetic field"
00015 << "\n-----------------------------------------------------------";
00016
00017 G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << G4endl;
00018
00019 std::ifstream file;
00020 file.open(filename);
00021 G4String exceptionString = "Unable to load field map file: " + (G4String)filename;
00022 if(!file) G4Exception(exceptionString.c_str(), "-1", FatalException, "");
00023
00024
00025 char buffer[256];
00026 file.getline(buffer,256);
00027
00028
00029 file >> nx >> ny >> nz;
00030
00031 G4cout << " [ Number of values x,y,z: "
00032 << nx << " " << ny << " " << nz << " ] "
00033 << G4endl;
00034
00035
00036 xField.resize( nx );
00037 yField.resize( nx );
00038 zField.resize( nx );
00039 int ix, iy, iz;
00040 for (ix=0; ix<nx; ix++) {
00041 xField[ix].resize(ny);
00042 yField[ix].resize(ny);
00043 zField[ix].resize(ny);
00044 for (iy=0; iy<ny; iy++) {
00045 xField[ix][iy].resize(nz);
00046 yField[ix][iy].resize(nz);
00047 zField[ix][iy].resize(nz);
00048 }
00049 }
00050
00051
00052
00053
00054 do {
00055 file.getline(buffer,256);
00056 } while ( buffer[1]!='0');
00057
00058
00059 double xval=0.0,yval=0.0,zval=0.0,bx,by,bz;
00060
00061 minx = miny = minz = 0.0;
00062 for (ix=0; ix<nx; ix++) {
00063 for (iy=0; iy<ny; iy++) {
00064 for (iz=0; iz<nz; iz++) {
00065 file >> xval >> yval >> zval >> bx >> by >> bz;
00066 G4cout << "Read: " << xval << " " << yval << " " << zval << " " << bx << " " << by << " " << bz << G4endl;
00067 if ( ix==0 && iy==0 && iz==0 ) {
00068 minx = xval * _lenUnit;
00069 miny = yval * _lenUnit;
00070 minz = zval * _lenUnit;
00071 }
00072 xField[ix][iy][iz] = bx * _fieldUnit;
00073 yField[ix][iy][iz] = by * _fieldUnit;
00074 zField[ix][iy][iz] = bz * _fieldUnit;
00075 }
00076 }
00077 }
00078 file.close();
00079
00080 maxx = xval * _lenUnit;
00081 maxy = yval * _lenUnit;
00082 maxz = zval * _lenUnit;
00083
00084 G4cout << "\n ---> ... done reading " << G4endl;
00085
00086
00087 G4cout << " ---> assumed the order: x, y, z, Bx, By, Bz "
00088 << "\n ---> Min values x,y,z: "
00089 << minx/CLHEP::cm << " " << miny/CLHEP::cm << " " << minz/CLHEP::cm << " cm "
00090 << "\n ---> Max values x,y,z: "
00091 << maxx/CLHEP::cm << " " << maxy/CLHEP::cm << " " << maxz/CLHEP::cm << " cm "
00092 << "\n ---> The field will be offset by " << zOffset/CLHEP::cm << " cm " << G4endl;
00093
00094
00095 if (maxx < minx) {std::swap(maxx,minx); invertX = true;}
00096 if (maxy < miny) {std::swap(maxy,miny); invertY = true;}
00097 if (maxz < minz) {std::swap(maxz,minz); invertZ = true;}
00098 G4cout << "\nAfter reordering if neccesary"
00099 << "\n ---> Min values x,y,z: "
00100 << minx/CLHEP::cm << " " << miny/CLHEP::cm << " " << minz/CLHEP::cm << " cm "
00101 << " \n ---> Max values x,y,z: "
00102 << maxx/CLHEP::cm << " " << maxy/CLHEP::cm << " " << maxz/CLHEP::cm << " cm ";
00103
00104 dx = maxx - minx;
00105 dy = maxy - miny;
00106 dz = maxz - minz;
00107 G4cout << "\n ---> Dif values x,y,z (range): "
00108 << dx/CLHEP::cm << " " << dy/CLHEP::cm << " " << dz/CLHEP::cm << " cm in z "
00109 << "\n-----------------------------------------------------------" << G4endl;
00110 }
00111
00112 void BDS3DMagField::GetFieldValue(const double point[4],
00113 double *Bfield ) const
00114 {
00115
00116 G4ThreeVector local;
00117
00118 local[0] = point[0] + translation[0];
00119 local[1] = point[1] + translation[1];
00120 local[2] = point[2] + translation[2] + fZoffset;
00121
00122 local *= Rotation();
00123
00124 #ifdef BDSDEBUG
00125 G4cout << "BDS3DMagField::GetFieldValue" << G4endl;
00126 G4cout << "point x = " << point[0]/CLHEP::cm << " cm" << G4endl;
00127 G4cout << "point y = " << point[1]/CLHEP::cm << " cm" << G4endl;
00128 G4cout << "point z = " << point[2]/CLHEP::cm << " cm" << G4endl;
00129 G4cout << "translation x = " << translation[0]/CLHEP::cm << " cm" << G4endl;
00130 G4cout << "translation y = " << translation[1]/CLHEP::cm << " cm" << G4endl;
00131 G4cout << "translation z = " << translation[2]/CLHEP::cm << " cm" << G4endl;
00132 G4cout << "fZOffset = " << fZoffset/CLHEP::cm << " cm" << G4endl;
00133 G4cout << "local x = " << local[0]/CLHEP::cm << " cm" << G4endl;
00134 G4cout << "local y = " << local[1]/CLHEP::cm << " cm" << G4endl;
00135 G4cout << "local z = " << local[2]/CLHEP::cm << " cm" << G4endl;
00136 #endif
00137
00138 double signy=1;
00139 double signz=1;
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 if ( local[0]>=minx && local[0]<=maxx &&
00163 local[1]>=miny && local[1]<=maxy &&
00164 local[2]>=minz && local[2]<=maxz ) {
00165
00166
00167
00168 double xfraction = (local[0] - minx) / dx;
00169 double yfraction = (local[1] - miny) / dy;
00170 double zfraction = (local[2] - minz) / dz;
00171
00172 if (invertX) { xfraction = 1 - xfraction;}
00173 if (invertY) { yfraction = 1 - yfraction;}
00174 if (invertZ) { zfraction = 1 - zfraction;}
00175
00176
00177
00178 double xdindex, ydindex, zdindex;
00179
00180
00181
00182 double xlocal = ( std::modf(xfraction*(nx-1), &xdindex));
00183 double ylocal = ( std::modf(yfraction*(ny-1), &ydindex));
00184 double zlocal = ( std::modf(zfraction*(nz-1), &zdindex));
00185
00186
00187
00188 int xindex = static_cast<int>(xdindex);
00189 int yindex = static_cast<int>(ydindex);
00190 int zindex = static_cast<int>(zdindex);
00191
00192
00193 #ifdef DEBUG_INTERPOLATING_FIELD
00194 G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << G4endl;
00195 G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << G4endl;
00196 double valx0z0, mulx0z0, valx1z0, mulx1z0;
00197 double valx0z1, mulx0z1, valx1z1, mulx1z1;
00198
00199
00200
00201
00202 #endif
00203
00204
00205 Bfield[0] =
00206 xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
00207 xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
00208 xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
00209 xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
00210 xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
00211 xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
00212 xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
00213 xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
00214 Bfield[1] = signy *(
00215 yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
00216 yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
00217 yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
00218 yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
00219 yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
00220 yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
00221 yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
00222 yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal );
00223 Bfield[2] = signz *(
00224 zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
00225 zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
00226 zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
00227 zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
00228 zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
00229 zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
00230 zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
00231 zField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal );
00232
00233 } else {
00234 Bfield[0] = 0.0;
00235 Bfield[1] = 0.0;
00236 Bfield[2] = 0.0;
00237 }
00238
00239 #ifdef BDSDEBUG
00240 G4cout << "Bfield x,y,z = " <<
00241 Bfield[0]/_fieldUnit << " " <<
00242 Bfield[1]/_fieldUnit << " " <<
00243 Bfield[2]/_fieldUnit <<
00244 G4endl;
00245 #endif
00246 }
00247
00248 void BDS3DMagField::Prepare(G4VPhysicalVolume *referenceVolume){
00249 const G4RotationMatrix* Rot= referenceVolume->GetFrameRotation();
00250 const G4ThreeVector Trans=referenceVolume->GetFrameTranslation();
00251 SetOriginRotation(*Rot);
00252 SetOriginTranslation(Trans);
00253 }