Hello Geant4 Gurus!!
I am using geant4.8.2.p01 currently and am implementing an electric field,
and later a magnetic field, to a specified local volume. I have defined
my own electric field class, all of the code is below for that class.
The main concern is point[4] which comes in at the GetFieldValue() function.
I have displayed point[0], point[1], point[2] and calculated the magnitude
of this vector and in some instances it is larger than my entire world volume,
which obviously makes no sense. I guess I am curious about the coordinate
system that point[4] is using... does it use the logicalVolume in which
the fieldManager of the field belongs to as the coordinate system, or is it
with respect to the global volume??
Also, when I place this field_vol into the global and it overlaps with other
volumes already there or new ones that are placed on, it does not give me a
'overlap' warning, and the protons (which is what I am using) go through
all volumes, overlapping or not, and it seems okay... is there something
special about a volume that is made of Vacuum and holds a field manager??
The main issue is what this point[4] is, because the EField is calculated
based on that number, so if I am using it incorrectly then there is a problem.
Thanks in advance for your help.
Cheers,
Mario
HEADER:
#ifndef ChargeShell_hh
#define ChargeShell_hh 1
#include "G4Types.hh"
#include "G4ThreeVector.hh"
#include "G4ElectricField.hh"
class ChargeShell : public G4ElectricField
{
public:
ChargeShell();
virtual ~ChargeShell() ;
virtual void GetFieldValue(const G4double point[4],
G4double *field) const;
G4ThreeVector CalcVolumeCharge(const G4double[4]) const;
G4ThreeVector CalcSurfaceCharge(const G4double[4]) const;
G4double GetLength() {return 10.0*ro;};
G4ThreeVector GetPos() {return srcToShellPos;};
private:
void ReadInputFile();
G4double ri;
G4double ro;
G4double Q;
G4ThreeVector srcToShellPos;
G4double fourPiEpsNot;
};
#endif
CC FILE:
#include "ChargeShell.hh"
extern int readFileVerbose;
ChargeShell::ChargeShell()
{
ReadInputFile();
fourPiEpsNot = 4.0*3.14159265359*8.85418781762e-12*farad/meter;
}
// ------------------------------------------------------------------------
ChargeShell::~ChargeShell()
{
}
void ChargeShell::GetFieldValue (const G4double point[4],G4double *fieldOut) const
{
if(ri > ro){
G4cout << "ri cannot be larger than ro" << G4endl;
return;
}
G4ThreeVector fieldVal = G4ThreeVector(0.0,0.0,0.0);
if(ro == ri){fieldVal = CalcSurfaceCharge(point);}
else{fieldVal = CalcVolumeCharge(point);}
G4cout << "fieldVal.x() = " << fieldVal.x()/(volt/m) << G4endl;
G4cout << "fieldVal.y() = " << fieldVal.y()/(volt/m) << G4endl;
G4cout << "fieldVal.z() = " << fieldVal.z()/(volt/m) << G4endl << G4endl;
fieldOut[0]= 0.0;
fieldOut[1]= 0.0;
fieldOut[2]= 0.0;
fieldOut[3]= fieldVal.x();
fieldOut[4]= fieldVal.y();
fieldOut[5]= fieldVal.z();
return;
}
G4ThreeVector ChargeShell::CalcSurfaceCharge(const G4double p[4]) const
{
G4double x = p[0];
G4double y = p[1];
G4double z = p[2];
G4double r = std::sqrt(x*x + y*y + z*z);
G4double Er;
G4cout << "r = " << r/mm << G4endl;
if(r <= ro){Er = 0.0*volt/m;G4cout << "----->>>>>Er = 0<<<<<-----" << G4endl;}
else{Er = Q/(fourPiEpsNot*r*r);}
G4cout << "Er = " << Er/(volt/m) << G4endl;
G4ThreeVector out = G4ThreeVector(Er*(x/r),Er*(y/r),Er*(z/r));
return out;
}
G4ThreeVector ChargeShell::CalcVolumeCharge(const G4double p[4]) const
{
G4double x = p[0];
G4double y = p[1];
G4double z = p[2];
G4double r = std::sqrt(x*x + y*y + z*z);
double riOverro = ri/ro;
double riOverr = ri/r;
G4double Er;
if(r >= 0.0 && r <= ri){Er = 0.0*volt/m;}
else if(r > ri && r < ro){
Er = (Q/(fourPiEpsNot*ro*ro))*(r/ro)*(1-riOverr*riOverr*riOverr)/(1-riOverro*riOverro*riOverro);
}
else{Er = Q/(fourPiEpsNot*r*r);}
G4ThreeVector out = G4ThreeVector(Er*(x/r),Er*(y/r),Er*(z/r));
return out;
}
void ChargeShell::ReadInputFile()
{
//read in parameters from text file
FILE *chargeShellInput;
char* inputFile = "C:\\g4work\\PRSim\\src\\Param_Files\\Param_ChargeShellField.txt";
chargeShellInput = fopen(inputFile,"rt");
if(chargeShellInput == NULL){
G4cout << "Param_ChargeShellField.txt did not come in correctly, check path" << G4endl;
G4cout << "chargeShellInput = " << chargeShellInput << G4endl;
G4cout << "errno = " << errno << G4endl;
return;
}
char readWholeLineBuff[500];
float tmpFloat;
char tmpString[50];
//read line of text for ri
fgets(readWholeLineBuff,500,chargeShellInput);
if(readWholeLineBuff != NULL){
fscanf(chargeShellInput,"%f\n",&tmpFloat);
ri = tmpFloat*um;
if(readFileVerbose > 0){
G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
G4cout << ri/um << G4endl;
}
}
else{G4cout << "Did not read in ri" << G4endl; return;}
//read line of text for ro
fgets(readWholeLineBuff,500,chargeShellInput);
if(readWholeLineBuff != NULL){
fscanf(chargeShellInput,"%f\n",&tmpFloat);
ro = tmpFloat*um;
if(readFileVerbose > 0){
G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
G4cout << ro/um << G4endl;
}
}
else{G4cout << "Did not read in ro" << G4endl; return;}
//read line of text for Q
fgets(readWholeLineBuff,500,chargeShellInput);
fgets(readWholeLineBuff,500,chargeShellInput);
if(readWholeLineBuff != NULL){
fscanf(chargeShellInput,"%f\n",&tmpFloat);
Q = tmpFloat*coulomb;
if(readFileVerbose > 0){
G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
G4cout << Q/coulomb << G4endl;
}
}
else{G4cout << "Did not read in Q" << G4endl; return;}
//read line of text for srcToShellPos
fgets(readWholeLineBuff,500,chargeShellInput);
if(readWholeLineBuff != NULL){
fscanf(chargeShellInput,"%f",&tmpFloat);
srcToShellPos.setX(tmpFloat*cm);
fscanf(chargeShellInput,"%f",&tmpFloat);
srcToShellPos.setY(tmpFloat*cm);
fscanf(chargeShellInput,"%f\n",&tmpFloat);
srcToShellPos.setZ(tmpFloat*cm);
if(readFileVerbose > 0){
G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
G4cout << srcToShellPos << G4endl;
}
}
else{G4cout << "Did not read in srcToShellPos" << G4endl; return;}
fclose(chargeShellInput);
return;
}
|
Local Electric Field: GetFieldValue() problem 
Feedback
]
icon
and there is no indentation. 
