| Message: Problem with magnetic field from a table | Not Logged In (login) |
|
Hi,
I have a 2d table with values of Bx and By (Bz is 0). It is quite large: 4900 points with spacing of 7.7mm. I use Geant 4.9.0.p01. I have a class:
---------------------VVV ------------------------
class WSMagneticField: public G4MagneticField
{
public:
WSMagneticField(G4ThreeVector); // The value of the field WSMagneticField(); // A zero field ~WSMagneticField();
void SetFieldValue(G4double fieldValue); void SetFieldValue(G4ThreeVector fieldVector);
void GetFieldValue(const double Point[3],
double *Bfield) const;
protected:
// Find the global Field Manager
G4FieldManager* GetGlobalFieldManager();
private:
static const int imesh=4900;
float x[imesh];
float y[imesh];
float Bx[imesh];
float By[imesh];
};
----------------------AAA---------------------
In the constructor I fill the tables x,y,Bx,By:
---------------------VVV----------------------
WSMagneticField::WSMagneticField()
{
// here open and read file
ifstream bmap;
string line;
const char *charline;
int ni1,ni2,idx=0;
float modB;
char ni3;
bmap.open("Bfield.data");
while(bmap) {
getline(bmap,line);
charline=line.c_str();
sscanf(charline," %d %d %f %f %f %f %f %c\n",
&ni1,&ni2,&x[idx],&y[idx],&Bx[idx],&By[idx],&modB,&ni3);
G4cout << "mag field check i = " << idx << " x = " << x[idx] << " Bx = " << Bx[idx] <<G4endl;
idx++;
}
GetGlobalFieldManager()->SetDetectorField(this); GetGlobalFieldManager()->CreateChordFinder(this); } --------------------AAA------------------------- and I have GetFieldValue function: --------------------VVV-------------------------
// argument is position in mm
void WSMagneticField::GetFieldValue(const double Point[3], double *Bfield) const
{
Bfield[0] = 0.;
Bfield[1] = 0.;
Bfield[2] = 0.;
G4float dmesh=7.93; // [mm]
G4double cPoint[3];
cPoint[0]=Point[0];
cPoint[1]=Point[1]-77.; // [mm] shift back the Quad to World ref.
cPoint[2]=Point[2];
G4double radius = pow(pow(cPoint[0],2)+pow(cPoint[1],2),0.5);
if(cPoint[0]>0 && radius < 240. && fabs(cPoint[2]) < 1400.) {
for(G4int i=0; i<3; i++) if(cPoint[i]<0) cPoint[i]*=(-1); // because of symmetry
for(G4int i=0; i<4900; i++) {
if(fabs(cPoint[0]-x[i])<dmesh/2 && fabs(cPoint[1]-y[i])<dmesh/2) {
Bfield[0]=Bx[i];
Bfield[1]=By[i];
break;
}
}
}
}
-------------------AAA-------------------------
And then, in DetectorConstruction I have:
-------------------VVV-------------------------
WSMagneticField* myField = new WSMagneticField();
G4FieldManager* fieldMgr
= G4TransportationManager::GetTransportationManager()
->GetFieldManager();
// I tried this and it did not work either // fieldMgr->SetDetectorField(myField); // fieldMgr->CreateChordFinder(myField); /* ... geometry declaration ... */
G4FieldManager *localFieldMgr = new G4FieldManager(myField); logicQ5->SetFieldManager(localFieldMgr,true); --------------------AAA------------------------- The program compiles, the file with table of values is read but never ever any event has finished being calculated. I've left it even for the whole night even if in normal conditions one event takes about 15 minutes. I would appreciate any suggestion leading to solution;-)
Mariusz Sapinski |
| Inline Depth: | Outline Depth: | Add message: |
|
to: |