365 lines
11 KiB
C++
365 lines
11 KiB
C++
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//*
|
|
// LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION
|
|
//
|
|
// ID :LEMuSRElectricField.cc , v 1.3
|
|
// AUTHOR: Taofiq PARAISO
|
|
// DATE : 2004-09-17 10:20
|
|
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//
|
|
//
|
|
// & &&&&&&&&&& &&&&&&& &&&&&&&&
|
|
// & & && && & &&
|
|
// & & & & & & &&
|
|
// & &&&&&&& & & &&&&&& &&&&&&&&
|
|
// & & & && & & &&
|
|
// & & && & & && && & &
|
|
// &&&&&&&&&& &&&&&&&&&& & &&&&& && &&&&&&& & &&
|
|
// &
|
|
// &
|
|
// &
|
|
// &
|
|
// Electric Field
|
|
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//
|
|
#include"LEMuSRElectricField.hh"
|
|
#include"G4UnitsTable.hh"
|
|
|
|
|
|
LEMuSRElectricField::LEMuSRElectricField(G4double fieldval,G4String Xfile,
|
|
G4String Yfile,
|
|
G4String Zfile,
|
|
G4String map_length_unit,
|
|
G4double Offset, G4double nbx, G4double nby, G4double nbz)
|
|
{
|
|
FieldVal=fieldval;
|
|
map_unit = map_length_unit;
|
|
G4cout<<G4endl
|
|
<< " ++++++ELECTRIC FIELD BEINg BUILT++++++"<<G4endl;
|
|
G4cout<< " +++++IMPORTANT++++++"<<G4endl;
|
|
G4cout<< "+ MAKE SURE UNITS IN FIELD MAP AND IN LEMuSRElectricField.cc ARE COHERENT +"<<G4endl<<G4endl;
|
|
|
|
// open files for reading
|
|
std::ifstream fx(Xfile);
|
|
std::ifstream fy(Yfile);
|
|
std::ifstream fz(Zfile);
|
|
|
|
// Ignore first blank line
|
|
char Xbuffer[256];
|
|
fx.getline(Xbuffer,256);
|
|
char Ybuffer[256];
|
|
fy.getline(Ybuffer,256);
|
|
char Zbuffer[256];
|
|
fz.getline(Zbuffer,256);
|
|
|
|
zOffset = Offset;
|
|
nx=nbx;
|
|
ny=nby;
|
|
nz=nbz;
|
|
|
|
|
|
int ix, iy, iz;
|
|
|
|
// Read in the data
|
|
double xval,yval,zval,bx,by,bz;
|
|
xval=0;
|
|
yval=0;
|
|
zval=0;
|
|
bx=0;
|
|
by=0;
|
|
bz=0;
|
|
ix=0;
|
|
iy=0;
|
|
iz=0;
|
|
G4int firstline;
|
|
firstline=0;
|
|
|
|
// ATTENTION :: FIELD MAP IS ASSUMED TO BEGIN WITH XMIN YMIN ZMIN AND END WITH XMAX YMAX ZMAX
|
|
// THIS IS A CASE FOR MAPS GENERATED BY FEMLAB
|
|
|
|
for (ix=0; ix<nx; ix++) {
|
|
for (iy=0; iy<ny; iy++) {
|
|
for (iz=0; iz<nz; iz++) {
|
|
|
|
fx >> xval >> yval >> zval >> bx;
|
|
fy >> xval >> yval >> zval >> by;
|
|
fz >> xval >> yval >> zval >> bz;
|
|
|
|
xField[ix][iy][iz] = bx*volt/meter;
|
|
yField[ix][iy][iz] = by*volt/meter;
|
|
zField[ix][iy][iz] = bz*volt/meter;
|
|
if (firstline==0)
|
|
{
|
|
minx =xval;
|
|
miny =yval;
|
|
minz =zval;
|
|
|
|
firstline=1;
|
|
}
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
fx.close(); fy.close(); fz.close();
|
|
|
|
maxx = xval;
|
|
maxy = yval;
|
|
maxz = zval;
|
|
|
|
G4cout<<G4endl
|
|
<< " (-: ++++++ELECTRIC FIELD INITIALIZED++++++ :-)";
|
|
G4cout<<"\n from files "<< Xfile << " etc. with following statistics : \n"
|
|
<<"\n last line: " << xval <<" "<< yval <<" "<< zval<<" "<<bx<<" "<<by<<" "<<bz<<" with field "<<fieldval<<"\n"
|
|
<<"\n boundaries: " << minx <<" "<< maxx <<" "<< miny<<" "<<maxy<<" "<<minz<<" "<<maxz<<" (in "<<map_unit<<")\n"
|
|
<<G4endl
|
|
<<G4endl;
|
|
|
|
|
|
//#ifdef DEBUG_INTERPOLATING_FIELD2
|
|
G4String output_name;
|
|
G4cout<<"Output Field Map file name?"<<G4endl;
|
|
G4cin>>output_name;
|
|
std::ofstream fout (output_name); // open destination file
|
|
fout.close ();
|
|
|
|
|
|
fout.open(output_name);
|
|
if (!fout.is_open()) exit(8);
|
|
|
|
for (ix=0; ix<nx; ix++) {
|
|
for (iy=0; iy<ny; iy++) {
|
|
for (iz=0; iz<nz; iz++) {
|
|
fout<< xField[ix][iy][iz]/volt*meter <<" " << yField[ix][iy][iz]/volt*meter <<" " << zField[ix][iy][iz]/volt*meter <<G4endl;
|
|
}
|
|
}
|
|
}
|
|
|
|
fout<< maxx<<" " << maxy<<" " <<maxz <<G4endl;
|
|
fout<< minx<<" " << miny<<" " <<minz <<G4endl;
|
|
|
|
fout.close();
|
|
|
|
G4cout<<"Field Map "<<output_name<<" genarated ___ press enter"<<G4endl;
|
|
|
|
//#endif
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
LEMuSRElectricField::LEMuSRElectricField(G4double fieldval,G4String file,
|
|
G4String map_length_unit,
|
|
G4double Offset,
|
|
G4double nbx,
|
|
G4double nby,
|
|
G4double nbz)
|
|
{
|
|
|
|
FieldVal= fieldval;
|
|
map_unit = map_length_unit;
|
|
G4cout<<G4endl
|
|
<< "++++++LEMUSR ELECTRIC FIELD INITIALIZATION+++++\n";
|
|
G4cout<< " +++++IMPORTANT++++++"<<G4endl;
|
|
G4cout<< "+ MAKE SURE UNITS IN FIELD MAP AND IN LEMuSRElectricField.cc ARE COHERENT +"<<G4endl;
|
|
G4cout<< "++++++ELECTRIC FIELD BEINg BUILT++++++"<<G4endl<<G4endl;
|
|
|
|
// open files for reading
|
|
std::ifstream fmap(file);
|
|
|
|
// Ignore first blank line
|
|
// char buffer[256];
|
|
// fmap.getline(buffer,256);
|
|
|
|
zOffset = Offset;
|
|
nx=nbx;
|
|
ny=nby;
|
|
nz=nbz;
|
|
|
|
int ix, iy, iz;
|
|
|
|
|
|
// Read in the data
|
|
double bx,by,bz;
|
|
for (ix=0; ix<nx; ix++) {
|
|
for (iy=0; iy<ny; iy++) {
|
|
for (iz=0; iz<nz; iz++) {
|
|
fmap >>bx >>by >>bz;
|
|
|
|
xField[ix][iy][iz] = bx*volt/meter;
|
|
yField[ix][iy][iz] = by*volt/meter;
|
|
zField[ix][iy][iz] = bz*volt/meter;
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
fmap>> maxx >> maxy >> maxz ;
|
|
fmap>> minx >> miny >> minz ;
|
|
G4cout<< maxx<<" " << maxy<<" " <<maxz <<G4endl;
|
|
G4cout<< minx<<" " << miny<<" " <<minz <<G4endl;
|
|
|
|
fmap.close();
|
|
|
|
G4cout<< "++++++ELECTRIC FIELD INITIALIZED++++++";
|
|
G4cout<<"\n from file "<< file <<G4endl<<G4endl;
|
|
|
|
}
|
|
|
|
|
|
LEMuSRElectricField::~LEMuSRElectricField()
|
|
{
|
|
delete [] xField ;
|
|
delete [] yField ;
|
|
delete [] zField ;
|
|
}
|
|
|
|
|
|
|
|
void LEMuSRElectricField::GetFieldValue(const G4double point[4],
|
|
G4double *Bfield ) const
|
|
{
|
|
// Values of x y z must be scaled according to the field map unit! important
|
|
// be carefull and use the debuging to verify it
|
|
|
|
G4double length_ratio, field_factor;
|
|
length_ratio=1;
|
|
field_factor=0;
|
|
|
|
if (map_unit=="m")
|
|
{
|
|
length_ratio=1000;
|
|
field_factor=1;
|
|
}
|
|
|
|
else if (map_unit=="dm")
|
|
{
|
|
length_ratio=100;
|
|
field_factor=10;
|
|
}
|
|
else if (map_unit=="cm")
|
|
{
|
|
length_ratio=10;
|
|
field_factor=100;
|
|
}
|
|
else if (map_unit=="mm")
|
|
{
|
|
length_ratio=1;
|
|
field_factor=1000;
|
|
}
|
|
else
|
|
{
|
|
G4cout<<"ERROR!!!! UNRECOGNIZED LENGTH UNIT FOR FIELD MAP ==>> zöro";
|
|
}
|
|
|
|
|
|
//! Scaling the distances acording to the unit.
|
|
G4double x = (point[0]/mm)/length_ratio;
|
|
G4double y = (point[1]/mm)/length_ratio;
|
|
G4double z = ((point[2] - zOffset)/mm)/length_ratio;
|
|
|
|
|
|
#ifdef DEBUG_INTERPOLATING_FIELD
|
|
|
|
G4cout<<"POSITION:: " <<point[0] <<" mm "<< G4BestUnit(point[1],"Length") <<" "<<G4BestUnit(point[2],"Length")<<G4endl;
|
|
|
|
G4cout<<"POSITION in the field map:: " <<x <<map_unit<<" "<< y <<map_unit<<" "<< z <<map_unit<< " "<<G4endl;
|
|
#endif
|
|
|
|
|
|
//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]
|
|
G4double xfraction = (x - minx) / (maxx-minx);
|
|
G4double yfraction = (y - miny) / (maxy-miny);
|
|
G4double zfraction = (z - minz) / (maxz-minz);//!!!!! PROBLEM MAY BE HERE
|
|
|
|
// Need addresses of these to pass to modf below.
|
|
// modf uses its second argument as an OUTPUT argument.
|
|
G4double xdindex, ydindex, zdindex;
|
|
|
|
// Position of the point within the cuboid defined by the
|
|
// nearest surrounding tabulated points
|
|
|
|
G4double xlocal = ( modf(xfraction*(nx-1), &xdindex));
|
|
G4double ylocal = ( modf(yfraction*(ny-1), &ydindex));
|
|
G4double zlocal = ( 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);
|
|
*/
|
|
G4int xindex = (G4int)(xdindex);
|
|
G4int yindex = (G4int)(ydindex);
|
|
G4int zindex = (G4int)(zdindex);
|
|
|
|
|
|
#ifdef DEBUG_INTERPOLATING_FIELD
|
|
G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << endl;
|
|
G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << endl;
|
|
#endif
|
|
|
|
//********************* FIELD INTERPOLATION***************************//
|
|
//*********** FROM GEANT4 EXTENDED EXAMPLE PURGING MAGNET ************//
|
|
|
|
// Full 3-dimensional version
|
|
|
|
// example: if map_unit= 'dm' multiply by 10 since field map generated by femlab is in volt/decimeter! in order to get the field in V/m.
|
|
|
|
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 {
|
|
// G4cout<<"RANGE ERROR=>field 0:: RANGE is " << minx <<" " << maxx <<" " << miny <<" " << maxy <<" " << minz <<" "<< maxz <<" ";
|
|
Bfield[0] = 0.0;
|
|
Bfield[1] = 0.0;
|
|
Bfield[2] = 0.0;
|
|
}
|
|
|
|
|
|
|
|
Bfield[0] = Bfield[0]*FieldVal*field_factor;
|
|
Bfield[1] = Bfield[1]*FieldVal*field_factor;
|
|
Bfield[2] = Bfield[2]*FieldVal*field_factor;
|
|
|
|
|
|
#ifdef DEBUG_INTERPOLATING_FIELD
|
|
G4cout<<"Field Value " << G4BestUnit(Bfield[0]*meter,"Electric potential") <<"/m " <<Bfield[1]/volt*meter <<" V/m "<< Bfield[2]/volt*meter <<" V/m " <<G4endl;
|
|
G4cout<<"FieldVal " <<FieldVal <<" coef " << field_factor<<" \n";
|
|
#endif
|
|
|
|
}
|
|
|
|
|
|
|