Dear Mariusz,
Geant4 uses the GetFieldValue method that you create many times during
every step. Looking at the method that you have written, I can see that
it is making a very large number of operations, in order to calculate
the value. So it is not surprising that it your program is very slow.
I have two suggestion to speed up your key GetFieldValue method:
1) Avoid large loops:
for(G4int i=0; i<4900; i++) {
if(fabs(cPoint[0]-x[i])<dmesh/2 ...
Use the functions for converting a floating point number to a
integer, in order to calculate the nearest 'i' integer to your point.
Some code along the lines of:
double ratio = (cPoint[0]/dmesh);
int i = round( ratio ); // Double check this - it
should round to nearest integer maybe it is std::round
if( fabs( cpoint[0]- x[i]) < dmesh/2 ) ...
etc
2) to replace the line
G4double radius = pow(pow(cPoint[0],2)+pow(cPoint[1],2),0.5);
with
G4double radius = std::sqrt(Point[0]*cPoint[0]+cPoint[1]*cPoint[1]);
since
- the power function is very costly compared to a multiplication, and
- square root is faster than pow( .. , 0.5)
The power method generally requires a logarithm and an exponential.
I expect that this will speed up your field method, and your
application, by a factor of 100 or more. I expect that your particles
will actually move.
Best regards,
John
Mariusz wrote:
> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
> In-Reply-To: <"/emfields/146"@geant4-hn.slac.stanford.edu>
> Subject: ...change this to be about your reply.
>
> 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
>
>
|
Re: Problem with magnetic field from a table 
Feedback
]
icon
and there is no indentation. 
