408 lines
14 KiB
C++
408 lines
14 KiB
C++
// Geant4 simulation for MuSR
|
|
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
|
|
// DATE : 2008-05
|
|
//
|
|
|
|
#include "lem4TabulatedElementField2D.hh"
|
|
//#include "lem4EventAction.hh" // cks: just to get the event nr. if needed.
|
|
#include "lem4Parameters.hh" /// Comment by TS - is this really needed?
|
|
///#include "lem4ErrorMessage.hh" -> PUT SOME LIMIT DURING WRONG UNIT CONVERSIONS
|
|
#include "G4UnitsTable.hh"
|
|
#include <string>
|
|
|
|
//lem4TabulatedElementField2D::lem4TabulatedElementField2D(const char* filename, G4double fieldValue, G4double lenUnit, G4double fieldNormalisation, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume),
|
|
// ffieldValue(fieldValue)
|
|
|
|
|
|
///G4double lenUnit, G4double fieldNormalisation,
|
|
|
|
lem4TabulatedElementField2D::lem4TabulatedElementField2D(const char* filename, const char fieldType, G4double fieldValue, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume),
|
|
fldType(fieldType), ffieldValue(fieldValue)
|
|
{
|
|
G4cout << "\n-----------------------------------------------------------";
|
|
|
|
// The DEFAULT user-defined field units for E and B (kilovolt/mm and tesla)
|
|
// NOTE: Should be the same as EMfieldUnit defined in DetectorConstruction.cc!
|
|
if (fldType == 'E') {G4cout << ("\n Electric field 2D"); fUnit = "kV/mm"; fieUnit= kilovolt/mm;}
|
|
else if (fldType == 'B') {G4cout << ("\n Magnetic field 2D"); fUnit = "T"; fieUnit = tesla;}
|
|
|
|
G4cout << "\n-----------------------------------------------------------" << G4endl;
|
|
G4cout << "\n ---> Reading the "<< fldType <<"-field map from " << filename << " ... " << G4endl;
|
|
|
|
// Open the file for reading
|
|
std::ifstream file( filename );
|
|
|
|
///char buffer[256]; //not needed with file peek
|
|
// Ignore first blank line
|
|
// file.getline(buffer,256);
|
|
|
|
// Read table dimensions
|
|
//G4String lUnit, fUnit;
|
|
//G4double cFact;
|
|
//fieldNormalisation = cFact;
|
|
//lenUnit = 1*cm,
|
|
//fieldNormalisation = 0.00001;
|
|
// 101 1 1001 cm tesla 0.00001
|
|
//G4double fieldNorm, lenNorm;
|
|
//G4cout<< lUnit << " XXXX " << fUnit << G4endl;
|
|
//G4cout<< "Tesla in G4 is "<< tesla << ", while kV/mm is " << kilovolt/mm << G4endl;
|
|
|
|
///file >> nr >> nDummy >> nz >> lUnit >> fUnit >> fieldNormalisation; // Note dodgy order
|
|
// Read the file header
|
|
file >> nr >> nz >> lUnit >> fieldNormalisation;
|
|
// Add manually the length unit and norm. factor in the field-map file if missing!
|
|
|
|
G4cout << " The grid consists of [" << nr << " x " << nz << "] r and z values" << G4endl;
|
|
|
|
// Get the field map LENGTH unit and its value in G4 native units.
|
|
lenUnit = G4UnitDefinition::GetValueOf(lUnit);
|
|
G4cout << " Field length unit = " << lUnit << ", in G4 native units = " << lenUnit << G4endl;
|
|
|
|
/*
|
|
// Get the EM field AMPLITUDE unit and its value in G4 native units.
|
|
// Separate the E case from B, since E presents complex units (e.g. kV/mm)
|
|
if (fldType == 'E') {
|
|
G4double volNorm, lenNorm;
|
|
std::string::size_type pos = fUnit.find("/");
|
|
std::string v_unit = fUnit.substr(0, pos); // Voltage unit
|
|
std::string l_unit = fUnit.substr(pos+1); // Length unit
|
|
|
|
volNorm = G4UnitDefinition::GetValueOf(v_unit);
|
|
lenNorm = G4UnitDefinition::GetValueOf(l_unit);
|
|
///G4cout<< v_unit << " " << l_unit << G4endl;
|
|
///G4cout<< volNorm << " " << lenNorm << G4endl;
|
|
|
|
fNorm = volNorm/lenNorm; // Electric field unit = [Voltage]/[Length]
|
|
G4cout << " Electric field unit = " << fUnit << ", in G4 native units = " << fNorm << G4endl;
|
|
}
|
|
|
|
else if (fldType == 'B') {
|
|
fNorm = G4UnitDefinition::GetValueOf(fUnit);
|
|
G4cout << " Magnetic field unit = " << fUnit << ", in G4 native units = " << fNorm << G4endl;
|
|
}
|
|
*/
|
|
G4cout << " Field map normalisation factor = " << fieldNormalisation << G4endl;
|
|
|
|
|
|
///fieldValue = fieldValue*fNorm/tesla; // Convert
|
|
|
|
///fieldValue = fieldValue/fNorm; // Convert
|
|
//double lenUnit = meter;
|
|
//double fieldUnit= tesla;
|
|
|
|
//lenNorm = G4UnitDefinition::GetValueOf(leng_unit);
|
|
//lenNorm = G4UnitDefinition::GetValueOf(lUnit);
|
|
//fieldNorm = G4UnitDefinition::GetValueOf(fNorm);
|
|
|
|
|
|
//G4double GetNewDoubleValue(const char*)
|
|
//fieldNorm = G4UIcmdWithADouble::SetNewDoubleValue(fNorm);
|
|
//fieldNorm = G4UIcmdWithADoubleAndUnit::GetNewUnitValue(fNorm); //GetNewDoubleRawValue(fNorm);
|
|
|
|
|
|
/*
|
|
|
|
G4cout << " Length Unit = " << leng_unit << G4endl; //lUnit
|
|
G4cout << " Field Unit = " << leng_unit << G4endl; //lUnit
|
|
|
|
G4cout << " Length Unit = " << leng_unit << G4endl; //lUnit
|
|
G4cout << " Conv. Unit = " << lenNorm << G4endl;
|
|
G4cout << " Field Unit = " << fNorm << G4endl;
|
|
G4cout << " Conv. Unit = " << fieldNorm << G4endl;
|
|
G4cout << " Norm. Fact = " << cFact << G4endl;
|
|
|
|
*/
|
|
/// FIELD NORMALISATION FACTOR INSIDE FIELD MAPS IS SUCH THAT (MAX. FIELD VALUE)*FNORMFACT = 1!
|
|
/// CHANGE SIGN TO REVERT FIELD DIRECTION! THE FIELD MAP HAS NO UNITS.
|
|
|
|
///G4cout << " Field set to "<< fieldValue/fNorm << " " << fUnit << G4endl; ///????
|
|
///G4cout << " Field set to "<< fieldValue*fNorm << " " << fUnit << G4endl;
|
|
///G4cout << " Field set to "<< fieldValue << " " << fUnit << G4endl;
|
|
|
|
|
|
// Set up storage space for the 2D 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.ignore(256, '\n');
|
|
} while (file.peek() == '%');
|
|
|
|
/// Attention: The old version works only ONCE, i.e. at a STOP signal, e.g. !
|
|
//do {
|
|
// file.getline(buffer, 256, '!');
|
|
//} while (buffer[0] != '!');
|
|
|
|
|
|
// Read in the data: [r, z, EMr, EMz]
|
|
double rval, zval, fr, fz;
|
|
for (ir=0; ir<nr; ir++) {
|
|
for (iz=0; iz<nz; iz++) {
|
|
file >> rval >> zval >> fr >> fz;
|
|
//G4cout << rval << " v " << zval << " v " << fr << " v " << fz << G4endl;
|
|
if ( ir==0 && iz==0 ) {
|
|
minimumr = rval * lenUnit;
|
|
minimumz = zval * lenUnit;
|
|
}
|
|
rField[ir][iz] = fr*fieldNormalisation;
|
|
zField[ir][iz] = fz*fieldNormalisation;
|
|
}
|
|
}
|
|
file.close();
|
|
|
|
maximumr = rval * lenUnit;
|
|
maximumz = zval * lenUnit;
|
|
|
|
G4String volumeName = logVolume->GetName().substr(4);
|
|
|
|
G4cout << " ... done reading " << G4endl;
|
|
G4cout << "\n ---> " << fldType << "-field in volume "<< volumeName
|
|
<< " set to: " << fieldValue/fieUnit << " " << fUnit << G4endl;
|
|
G4cout << "\n ---> Assumed order: r, z, "<< fldType <<"r, "<< fldType <<"z "
|
|
<< "\n Min values r, z: "
|
|
<< minimumr/cm << " " << minimumz/cm << " cm "
|
|
<< "\n Max values r, z: "
|
|
<< maximumr/cm << " " << maximumz/cm << " cm" << G4endl;
|
|
|
|
|
|
// Should really check that the limits are not the wrong way around.
|
|
if (maximumr < minimumr) {Invert("r");}
|
|
if (maximumz < minimumz) {Invert("z");}
|
|
|
|
if ((maximumr < minimumr) || (maximumz < minimumz)) {
|
|
G4cout << "\n ---> After reordering:"
|
|
<< "\n Min values r, z: "
|
|
<< minimumr/cm << " " << minimumz/cm << " cm "
|
|
<< "\n Max values r, z: "
|
|
<< maximumr/cm << " " << maximumz/cm << " cm " << G4endl;
|
|
}
|
|
|
|
dr = maximumr - minimumr;
|
|
dz = maximumz - minimumz;
|
|
|
|
G4cout << " Range of values: "
|
|
<< dr/cm << " cm (in r) and " << dz/cm << " cm (in z)"
|
|
<< "\n-----------------------------------------------------------\n" << G4endl;
|
|
|
|
}
|
|
|
|
|
|
void lem4TabulatedElementField2D::addFieldValue(const G4double point[4],
|
|
G4double *field ) const
|
|
{
|
|
|
|
///G4cout << "addFieldValue is being called" << G4endl;
|
|
G4double EMf[3]; // EM field value obtained from the field map
|
|
///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 = local.z(); //Modified by TS
|
|
/// double z = fabs(local.z());
|
|
/// double z_sign =(local.z()>0) ? 1.:-1.;
|
|
|
|
// G4cout<<"Global points= "<<point[0]<<", "<<point[1]<<", "<<point[2]<<", Local point= "<<x<<", "<<y<<", "<<z<<G4endl;
|
|
|
|
// Check that the point is within the defined region
|
|
if ( r<maximumr && z<maximumz ) {
|
|
// 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 - minimumr) / dr;
|
|
double zfraction = (z - minimumz) / 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 rindex 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 EMfield_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 ;
|
|
EMf[0] = (r>0) ? EMfield_R * (local.x() /r) : 0.;
|
|
EMf[1] = (r>0) ? EMfield_R * (local.y() /r) : 0.;
|
|
EMf[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 ;
|
|
|
|
EMf[0] *= ffieldValue; /// * z_sign; ///Modified by TS: Not need for FULL field map
|
|
EMf[1] *= ffieldValue; /// * z_sign; ///Modified by TS: Not need for FULL field map
|
|
EMf[2] *= ffieldValue;
|
|
|
|
///G4cout << "The local coord. is: x_l=" << rlocal << ", z_l=" << zlocal << G4endl;
|
|
///G4cout << "The coord. is: x_loc=" << local.x() << ", y_loc=" << local.y() << ", z_loc=" << local.z() << G4endl;
|
|
|
|
G4ThreeVector finalField(EMf[0],EMf[1],EMf[2]);
|
|
finalField = global2local.Inverse().TransformAxis(finalField);
|
|
|
|
if (fldType == 'B') {
|
|
field[0] += finalField.x();
|
|
field[1] += finalField.y();
|
|
field[2] += finalField.z();
|
|
}
|
|
else if (fldType == 'E') {
|
|
field[3] += finalField.x();
|
|
field[4] += finalField.y();
|
|
field[5] += finalField.z();
|
|
}
|
|
|
|
}
|
|
// G4cout<<"Kamil: Field: ("<<field[0]/tesla<<","<<field[1]/tesla<<","<<field[2]/tesla<<")"<<G4endl;
|
|
|
|
}
|
|
|
|
|
|
|
|
G4double lem4TabulatedElementField2D::GetNominalFieldValue() {
|
|
return ffieldValue;
|
|
}
|
|
|
|
void lem4TabulatedElementField2D::SetNominalFieldValue(G4double newFieldValue) {
|
|
// Rescale the magnetic field for a new value of the magnetic field
|
|
ffieldValue=newFieldValue;
|
|
|
|
///G4cout<<"lem4TabulatedElementField2D.cc: ffieldValue changed to "<< ffieldValue << fUnit << G4endl;
|
|
G4cout<<"lem4TabulatedElementField2D.cc: ffieldValue changed to "<< ffieldValue/fieUnit << fUnit << G4endl;
|
|
///G4cout<<"lem4TabulatedElementField2D.cc: ffieldValue changed to "<< ffieldValue/tesla << " T"<< G4endl;
|
|
}
|
|
|
|
void lem4TabulatedElementField2D::Invert(const char* indexToInvert) {
|
|
// This function inverts the indices of the field table for a given axis (r or z).
|
|
// It should be called in the case when the r or z coordinate in the initial
|
|
// field table is ordered in the 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 lem4TabulatedElementField2D::Invert() function works properly!"<<G4endl;
|
|
G4cout<<"It has not been tested yet!"<<G4endl;
|
|
|
|
if (strcmp(indexToInvert,"r")==0) {invertR=true; std::swap(maximumr,minimumr);}
|
|
if (strcmp(indexToInvert,"z")==0) {invertZ=true; std::swap(maximumz,minimumz);}
|
|
|
|
for (int ir=0; ir<nr; ir++) {
|
|
for (int iz=0; iz<nz; iz++) {
|
|
if (invertR) {
|
|
rField[ir][iz] = rFieldTemp[nr-1-ir][iz];
|
|
zField[ir][iz] = zFieldTemp[nr-1-ir][iz];
|
|
}
|
|
else if(invertZ) {
|
|
rField[ir][iz] = rFieldTemp[ir][nz-1-iz];
|
|
zField[ir][iz] = zFieldTemp[ir][nz-1-iz];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
//#include "G4UIcmdWithADouble.hh"
|
|
//#include "G4UIcmdWithADoubleAndUnit.hh"
|
|
|
|
//#include <iostream>
|
|
//using namespace std;
|
|
|
|
///G4double G4UIcmdWithADouble::GetNewDoubleValue(const char* paramString)
|
|
|
|
|
|
///**********************///
|
|
//char str[] ="- This, a sample string.";
|
|
///char str[] ="kilovolt/mm";
|
|
///char * pch;
|
|
///printf ("Splitting string \"%s\" into tokens:\n",str);
|
|
//pch = strtok (str," ,.-");
|
|
///pch = strtok (str,"*/");
|
|
///while (pch != NULL)
|
|
///{
|
|
///printf ("%s\n",pch);
|
|
//pch = strtok (NULL, " ,.-");
|
|
///pch = strtok (NULL,"*/");
|
|
///}
|
|
|
|
///**********************///
|
|
|
|
//char str[] ="kilovolt/mm";
|
|
//char * pch;
|
|
//printf ("Splitting string \"%s\" into tokens:\n",str);
|
|
//pch = strtok (str," ,.-");
|
|
//pch = strtok (str,"/");
|
|
//for (int i=0; i<2; x++) {
|
|
//while (pch != NULL)
|
|
//{
|
|
// printf ("%s\n",pch);
|
|
// pch[]
|
|
//}
|
|
//pch = strtok (NULL,"*/");
|
|
//}
|
|
|
|
//int x;
|
|
//for(x=1; x <= 100; x++) {
|
|
//printf("%d ", x);
|
|
//}
|
|
|
|
///**********************///
|
|
// Only if field is E!
|
|
///char str[] = "This is a sample string kilovolt/mm";
|
|
///char * pch;
|
|
///printf ("Looking for the '/' character in \"%s\"...\n",str);
|
|
///pch = strchr(str,'/');
|
|
///while (pch!=NULL)
|
|
///{
|
|
///printf ("found at %d\n",pch-str+1);
|
|
///pch=strchr(pch+1,'s');
|
|
///}
|
|
|
|
///**********************///
|
|
|