347 lines
13 KiB
C++
347 lines
13 KiB
C++
// Geant4 simulation for MuSR
|
|
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
|
|
// DATE : 2008-05
|
|
//
|
|
|
|
#include "lem4Axial2DElField.hh"
|
|
#include "G4UnitsTable.hh"
|
|
//#include "lem4Parameters.hh"
|
|
|
|
lem4Axial2DElField::lem4Axial2DElField(const char* filename, G4double fieldValue, G4double lenUnit, G4double fieldNormalisation, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume),
|
|
ffieldValue(fieldValue)
|
|
|
|
//lem4Axial2DElField::lem4Axial2DElField(const G4String filename, double fieldValue, double lenUnit, double fieldNormalisation, double offset) : ffieldValue(fieldValue/1000.0), zOffset(offset), invertR(false), invertZ(false)
|
|
{
|
|
// double lenUnit = decimeter; = 100 // Base unit mm
|
|
// double fieldUnit = kV; = 0.001 // Base unit MV
|
|
// double fieldNorm = 0.00001 = 1E-5 // = 0.001/100
|
|
G4cout << "\n-----------------------------------------------------------"
|
|
<< "\n Electric field"
|
|
<< "\n-----------------------------------------------------------"
|
|
<< G4endl;
|
|
// G4cout << " kilovolt = "<< kilovolt << G4endl; // = 1e-3 MV - default G4 unit
|
|
///G4cout << "\n Potential set to "<< fieldValue/1000.0/kilovolt << " kV"<< G4endl;
|
|
G4cout << "\n Potential set to "<< fieldValue/kilovolt << " kV"<< G4endl;
|
|
G4cout << "\n ---> " "Reading 2D E-field map from " << filename << " ... " << G4endl;
|
|
std::ifstream file(filename); // Open the file for reading
|
|
|
|
// Start reading file content
|
|
char buffer[256];
|
|
|
|
// Read table dimensions first
|
|
file >> nr >> nz; // r - radial, z - longitudinal coordinate
|
|
|
|
G4cout << " 2D field table dimensions (r-z): ["
|
|
<< nr << ", " << nz << "] "
|
|
<< G4endl;
|
|
|
|
// Set up storage space for table
|
|
rField.resize(nr);
|
|
zField.resize(nr);
|
|
int ir, iz;
|
|
|
|
for (ir = 0; ir < nr; ir++) {
|
|
rField[ir].resize(nz);
|
|
zField[ir].resize(nz);
|
|
}
|
|
|
|
// Ignore header information. All lines whose first character
|
|
// is '%' are considered to be part of the header.
|
|
do {
|
|
file.getline(buffer, 256);
|
|
} while (buffer[0] != '%');
|
|
|
|
// Read in the data: [r, z, Er, Ez]
|
|
double rval, zval, er, ez;
|
|
for (ir = 0; ir < nr; ir++) {
|
|
for (iz = 0; iz < nz; iz++) {
|
|
file >> rval >> zval >> er >> ez;
|
|
if (ir == 0 && iz == 0) {
|
|
minr = rval * lenUnit;
|
|
minz = zval * lenUnit;
|
|
}
|
|
rField[ir][iz] = er*fieldNormalisation;
|
|
zField[ir][iz] = ez*fieldNormalisation;
|
|
}
|
|
}
|
|
G4cout << " Normalisation coeff.: " << fieldNormalisation << G4endl;
|
|
file.close();
|
|
|
|
maxr = rval * lenUnit;
|
|
maxz = zval * lenUnit;
|
|
|
|
G4cout << "\n ---> ... done reading (assumed order: r, z, Er, Ez)." << G4endl;
|
|
G4cout << "\n ---> Min values r, z: "
|
|
<< minr/cm << " " << minz/cm << " cm "
|
|
<< "\n ---> Max values r, z: "
|
|
<< maxr/cm << " " << maxz/cm << " cm " << G4endl;
|
|
|
|
// Should really check that the limits are not the wrong way around.
|
|
if (maxr < minr) {Invert("x");} ///{std::swap(maxr, minr); invertR = true;}
|
|
if (maxz < minz) {Invert("z");} ///{std::swap(maxz, minz); invertZ = true;}
|
|
if (maxr < minr || maxz < minz){ ///(invertR == true || invertZ == true){
|
|
G4cout << "\n After reordering:"
|
|
<< "\n ---> Min values r, z: "
|
|
<< minr/cm << " " << minz/cm << " cm "
|
|
<< "\n ---> Max values r, z: "
|
|
<< maxr/cm << " " << maxz/cm << " cm " << G4endl;
|
|
}
|
|
|
|
dr = maxr - minr;
|
|
dz = maxz - minz;
|
|
|
|
G4cout << " ---> Range of values: "
|
|
<< dr/cm << " cm (in r) and " << dz/cm << " cm (in z)."
|
|
<< "\n-----------------------------------------------------------\n" << G4endl;
|
|
}
|
|
|
|
|
|
///lem4Axial2DElField::~lem4Axial2DElField() /// Made virtual!
|
|
///{
|
|
/// delete [] &rField ;
|
|
/// delete [] &zField ;
|
|
///}
|
|
|
|
|
|
void lem4Axial2DElField::addFieldValue(const G4double point[4],
|
|
G4double *field) const
|
|
{
|
|
G4double B[3]; // Field value obtained from the field table
|
|
///TS
|
|
///G4cout<<"TTfileTT Global coord. ("<<point[0]/mm<<", "<<point[1]/mm<<", "<<point[2]/mm<<")."<<G4endl;
|
|
G4ThreeVector global(point[0],point[1],point[2]);
|
|
G4ThreeVector local;
|
|
|
|
local = global2local.TransformPoint(global);
|
|
///TS
|
|
///G4cout<<"TTfileTT Local coord. ("<<local[0]/mm<<", "<<local[1]/mm<<", "<<local[2]/mm<<")."<<G4endl;
|
|
///G4cout<<"The LOGICAL volume is named "<<lvolume->GetName()<<G4endl;
|
|
|
|
double r = sqrt(local.x()*local.x()+local.y()*local.y());
|
|
double z = fabs(local.z());
|
|
double z_sign = (local.z()>0) ? 1.:-1.;
|
|
|
|
// Check that the point is within the defined region
|
|
if ( r < maxr && z < maxz ) {
|
|
// Relative point position within region, normalized to [0,1].
|
|
double rfraction = (r - minr) / dr;
|
|
double zfraction = (z - minz) / dz;
|
|
|
|
// Need addresses of these to pass to modf below.
|
|
// modf uses its second argument as an OUTPUT argument.
|
|
double rdindex, zdindex;
|
|
|
|
// Position of the point within the cuboid defined by the
|
|
// nearest surrounding tabulated points
|
|
double rlocal = ( modf(rfraction*(nr-1), &rdindex));
|
|
double 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 rindex = static_cast<int>(rdindex);
|
|
int zindex = static_cast<int>(zdindex);
|
|
|
|
//cks The following check is necessary - even though xindex and zindex should never be out of range,
|
|
// it may happen (due to some rounding error ?). It is better to leave the check here.
|
|
if ((rindex<0)||(rindex>(nr-2))) {
|
|
std::cout<<"SERIOUS PROBLEM: rindex out of range! rindex="<<rindex<<" r="<<r<<" rfraction="<<rfraction<<std::endl;
|
|
if (rindex<0) rindex=0;
|
|
else rindex=nr-2;
|
|
}
|
|
if ((zindex<0)||(zindex>(nz-2))) {
|
|
std::cout<<"SERIOUS PROBLEM: zindex out of range! zindex="<<zindex<<" z="<<z<<" zfraction="<<zfraction<<std::endl;
|
|
if (zindex<0) zindex=0;
|
|
else zindex=nz-2;
|
|
}
|
|
|
|
// G4cout<<"xField["<<xindex<<"]["<<zindex<<"]="<<xField[xindex ][zindex ]<<G4endl;
|
|
// G4cout<<"zField["<<xindex<<"]["<<zindex<<"]="<<zField[xindex ][zindex ]<<G4endl;
|
|
|
|
// Interpolate between the neighbouring points
|
|
double Bfield_R =
|
|
rField[rindex ][zindex ] * (1-rlocal) * (1-zlocal) +
|
|
rField[rindex ][zindex+1] * (1-rlocal) * zlocal +
|
|
rField[rindex+1][zindex ] * rlocal * (1-zlocal) +
|
|
rField[rindex+1][zindex+1] * rlocal * zlocal ;
|
|
B[0] = (r>0) ? Bfield_R * (local.x() /r) : 0.;
|
|
B[1] = (r>0) ? Bfield_R * (local.y() /r) : 0.;
|
|
B[2] =
|
|
zField[rindex ][zindex ] * (1-rlocal) * (1-zlocal) +
|
|
zField[rindex ][zindex+1] * (1-rlocal) * zlocal +
|
|
zField[rindex+1][zindex ] * rlocal * (1-zlocal) +
|
|
zField[rindex+1][zindex+1] * rlocal * zlocal ;
|
|
|
|
B[0] *= ffieldValue * z_sign;
|
|
B[1] *= ffieldValue * z_sign;
|
|
B[2] *= ffieldValue;
|
|
|
|
G4ThreeVector finalField(B[0],B[1],B[2]);
|
|
finalField = global2local.Inverse().TransformAxis(finalField);
|
|
|
|
field[0] += finalField.x();
|
|
field[1] += finalField.y();
|
|
field[2] += finalField.z();
|
|
}
|
|
// G4cout<<"Kamil: Field: ("<<field[0]/tesla<<","<<field[1]/tesla<<","<<field[2]/tesla<<")"<<G4endl;
|
|
|
|
}
|
|
|
|
|
|
|
|
/* OLD IMPLEMENTATION - CHANGED WITH GLOBALFIELD! TS
|
|
void lem4Axial2DElField::GetFieldValue(const double point[4], double *Efield) const
|
|
{
|
|
// musrEventAction* myEventAction= musrEventAction::GetInstance();
|
|
// G4int evNr=myEventAction->GetLatestEventNr();
|
|
// G4int evNrKriz=6795; //457; //14250
|
|
// if (evNr==evNrKriz) {
|
|
// std::cout<<"evNr="<<evNr<<std::endl;
|
|
// printf ("for point= %f %f %f B= %10.10f %10.10f %10.10f \n",
|
|
// point[0],point[1],point[2],Bfield[0]/tesla,Bfield[1]/tesla,Bfield[2]/tesla);
|
|
// }
|
|
Efield[0]=0.; Efield[1]=0.; Efield[2]=0.;
|
|
|
|
double r = sqrt(point[0]*point[0]+point[1]*point[1]);
|
|
//Apply a fixed offset (in mm) to the z - coordinate.
|
|
double rel_z = (point[2] - zOffset)/mm;
|
|
double z = fabs(rel_z);
|
|
double z_sign = (rel_z>0) ? 1.:-1.;
|
|
//if (evNr==evNrKriz) std::cout<<"r ="<<r<<" maxr="<<maxr<<std::endl;
|
|
|
|
#ifdef DEBUG_INTERPOLATING_FIELD
|
|
double dst = sqrt(point[0]*point[0] + point[1]*point[1]);
|
|
G4cout<<"POSITION: "<< G4BestUnit(point[0], "Length")<<" "<< G4BestUnit(point[1], "Length") <<" "<< G4BestUnit(dst, "Length") <<" "<< G4BestUnit(point[2], "Length") <<G4endl;
|
|
G4cout<<"POSITION in field map: "<< G4BestUnit(r, "Length") <<" "<< G4BestUnit(z, "Length")<<G4endl;
|
|
#endif
|
|
|
|
|
|
// Check that the point is within the defined region
|
|
if ( r<maxr && z<maxz ) {
|
|
// if (evNr>evNrKriz) std::cout<<"bol som tu"<<std::endl;
|
|
|
|
// Position of given point within region, normalized to the range
|
|
// [0,1]
|
|
double rfraction = (r - minr) / dr;
|
|
double zfraction = (z - minz) / dz;
|
|
|
|
if (invertR) { rfraction = 1 - rfraction;}
|
|
if (invertZ) { zfraction = 1 - zfraction;}
|
|
// Need addresses of these to pass to modf below.
|
|
// modf uses its second argument as an OUTPUT argument.
|
|
double rdindex, zdindex;
|
|
|
|
// Position of the point within the cuboid defined by the
|
|
// nearest surrounding tabulated points
|
|
double rlocal = ( modf(rfraction*(nr-1), &rdindex));
|
|
double 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 rindex = static_cast<int>(rdindex);
|
|
int zindex = static_cast<int>(zdindex);
|
|
|
|
//cks The following check is necessary - even though xindex and zindex should never be out of range,
|
|
// it may happen (due to some rounding error ?). It is better to leave the check here.
|
|
if ((rindex<0)||(rindex>(nr-2))) {
|
|
std::cout<<"SERIOUS PROBLEM: rindex out of range! rindex="<<rindex<<" r="<<r<<" rfraction="<<rfraction<<std::endl;
|
|
if (rindex<0) rindex=0;
|
|
else rindex=nr-2;
|
|
}
|
|
if ((zindex<0)||(zindex>(nz-2))) {
|
|
std::cout<<"SERIOUS PROBLEM: zindex out of range! zindex="<<zindex<<" z="<<z<<" zfraction="<<zfraction<<std::endl;
|
|
if (zindex<0) zindex=0;
|
|
else zindex=nz-2;
|
|
}
|
|
|
|
#ifdef DEBUG_INTERPOLATING_FIELD
|
|
G4cout << "Local r, z: " << rlocal << " " << zlocal << G4endl;
|
|
G4cout << "Index r, z: " << rindex << " " << zindex << G4endl;
|
|
#endif
|
|
|
|
// Interpolate between the neighbouring points
|
|
double Efield_R =
|
|
rField[rindex ][zindex ] * (1-rlocal) * (1-zlocal) +
|
|
rField[rindex ][zindex+1] * (1-rlocal) * zlocal +
|
|
rField[rindex+1][zindex ] * rlocal * (1-zlocal) +
|
|
rField[rindex+1][zindex+1] * rlocal * zlocal ;
|
|
Efield[0] = (r>0) ? Efield_R * (point[0]/r) : 0.;
|
|
Efield[1] = (r>0) ? Efield_R * (point[1]/r) : 0.;
|
|
Efield[2] =
|
|
zField[rindex ][zindex ] * (1-rlocal) * (1-zlocal) +
|
|
zField[rindex ][zindex+1] * (1-rlocal) * zlocal +
|
|
zField[rindex+1][zindex ] * rlocal * (1-zlocal) +
|
|
zField[rindex+1][zindex+1] * rlocal * zlocal ;
|
|
// ATTENTION if dealing with MAGNETIC FIELDS reverse x, y and z signs!!
|
|
Efield[0] = Efield[0] * ffieldValue; // * z_sign;
|
|
Efield[1] = Efield[1] * ffieldValue; // * z_sign;
|
|
Efield[2] = Efield[2] * ffieldValue * z_sign;
|
|
}
|
|
|
|
|
|
|
|
// Set some small field if field is almost zero (to avoid internal problems of Geant).
|
|
if (sqrt(Efield[0]*Efield[0]+Efield[1]*Efield[1]+Efield[2]*Efield[2])<1.0E-12*volt) {
|
|
// if (evNr>evNrKriz) std::cout<<"bol som tuna"<<std::endl;
|
|
// Bfield[0] = 0.0;
|
|
// Bfield[1] = 0.0;
|
|
// Bfield[2] = Bfield[2] + 0.00001*tesla;
|
|
#ifdef DEBUG_INTERPOLATING_FIELD
|
|
G4cout<<"Zero field case"<<G4endl;
|
|
#endif
|
|
Efield[2] = 1.0E-12*volt;
|
|
}
|
|
// Print the field (for debugging)
|
|
// if ((point[2]>-1*cm)&&(point[2]<1*cm)) {
|
|
// printf ("for point= %f %f %f B= %10.10f %10.10f %10.10f \n",
|
|
// point[0],point[1],point[2],Bfield[0]/tesla,Bfield[1]/tesla,Bfield[2]/tesla);
|
|
// }
|
|
// if (evNr>evNrKriz) std::cout<<"this is the end"<<std::endl;
|
|
|
|
#ifdef DEBUG_INTERPOLATING_FIELD
|
|
//G4cout << " Volt/meter factor = "<< volt/meter << G4endl; // = 1e-9
|
|
G4cout<<"Field Value " << G4BestUnit(Efield[0]*meter,"Electric potential") <<"/m " <<Efield[1]/volt*meter <<" V/m "<< Efield[2]/volt*meter <<" V/m " <<G4endl;
|
|
G4cout<<"FieldVal " << ffieldValue/kilovolt << G4endl;
|
|
#endif
|
|
|
|
}
|
|
|
|
*/
|
|
|
|
G4double lem4Axial2DElField::GetNominalFieldValue() {
|
|
return ffieldValue;
|
|
}
|
|
|
|
void lem4Axial2DElField::SetNominalFieldValue(G4double newFieldValue) {
|
|
// Rescale the EM field to account for the new applied field value
|
|
ffieldValue=newFieldValue;
|
|
G4cout<<"lem4Axial2DElField.cc: ffieldValue changed to = "<< ffieldValue/kilovolt<<" kV."<<G4endl;
|
|
}
|
|
|
|
void lem4Axial2DElField::Invert(const char* indexToInvert) {
|
|
// This function inverts the field table indices along a given axis (r or z).
|
|
// It should be called whenever x or z coordinates in the initial field table
|
|
// are ordered in a decreasing order.
|
|
std::vector< std::vector< double > > rFieldTemp(rField);
|
|
std::vector< std::vector< double > > zFieldTemp(zField);
|
|
G4bool invertR=false;
|
|
G4bool invertZ=false;
|
|
|
|
G4cout<<"Check that the lem4Axial2DElField::Invert() function works properly!"<<G4endl;
|
|
G4cout<<"It has not been tested yet!"<<G4endl;
|
|
|
|
if (strcmp(indexToInvert,"r")==0) {invertR=true; std::swap(maxr, minr);}
|
|
if (strcmp(indexToInvert,"z")==0) {invertZ=true; std::swap(maxz, minz);}
|
|
|
|
for (int ir=0; ir<nr; ir++) {
|
|
for (int iz=0; iz<nz; iz++) {
|
|
if (invertR) {
|
|
rField[ir][iz] = rFieldTemp[nr-1-ir][iz];
|
|
}
|
|
else if(invertZ) {
|
|
rField[ir][iz] = rFieldTemp[ir][nz-1-iz];
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|