Message: Re: The Strangest Problem I've ever met Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

Question Re: The Strangest Problem I've ever met 

Forum: Fields: Magnetic and Otherwise
Re: Question The Strangest Problem I've ever met (Zukai Wang)
Re: Note Re: The Strangest Problem I've ever met (Tom Roberts)
Date: 20 Feb, 2009
From: Zukai Wang <Zukai Wang>

On Thu, 19 Feb 2009 15:24:38 GMT, Tom Roberts wrote:
> You most definitely do not want to open and read a file during EVERY
> call to GetFieldValue() -- this will slow your simulation down by a
> truly enormous factor! You consider "only 3 or 4 min" to be "quickly" --
> I routinely track hundreds of particles per second for many meters
> through magnetic field maps. Even a single sprintf() or sscanf() during
> tracking slows it down unacceptably.
> 
> Arrange the code so you open and read the files ONCE during
> initialization and keep the results around to use in GetFieldValue(). It
> is generally better (more general, easier to modify) to use a single
> file with a field map, rather than multiple files. A good design has the
> field map file describe its own geometrical structure (boundaries of the
> map, grid within the map, units, etc.) and the code uses this rather
> than hard-coded values for the geometry. I strongly recommend an ASCII
> format, so you can edit it easily, and can write conversion programs
> from the output of an EM modeling program (superfish, tosca, ...) to the
> format your simulation reads. Field maps of several tens of megabytes
> are not a problem on modern computers.
> 
> Whenever a track reaches the edge of the world, the track is killed.
> Similarly, when a track stops (kinetic energy goes to 0 or below a limit
> set in the program) it is killed unless there is an applicable physics
> process that makes it get captured or decay (which also kills the
> track). Just think about it -- how could a track possibly be tracked
> outside the known world? Or tracked forever while stopped?
> 
> Normally, when all tracks have been processed (i.e. tracked until they
> are killed), the run is ended and the program exits. But you must do
> this in your main program, as in the examples.
> 
Thank you so much for the direction. I have modified my code this way:
------------------------------------------------------------------------
ExN02MagneticField::ExN02MagneticField()
{
  fstream in;
  G4int n=0;
 in.open("Bmap.dat");

  while (in.good()&&!in.eof()) {
  in>>  x[n]>>y[n]>>z[n]>>Bx[n]>>By[n]>>Bz[n];

 G4cout<<  x[n]<<  y[n]<<  z[n]<<  Bx[n]<<  By[n]<<  Bz[n]<< G4endl;
 n++;}
  in.close();
  GetGlobalFieldManager()->SetDetectorField(this);
  GetGlobalFieldManager()->CreateChordFinder(this);

 }
ExN02MagneticField::~ExN02MagneticField()
{;}

void ExN02MagneticField::GetFieldValue(const double Point[3],double *Bfield) const
{
     Bfield[0]=0.;
     Bfield[1]=0.;
     Bfield[2]=0.;
    //G4int number=0;
// if (fabs(Point[0])<75.*mm && fabs(Point[1])<405.*mm && fabs(Point[2])<705.*mm)
  // Bfield[0]=-6450.*gauss;
    for (G4int i=0;i<23288;i++) 
{
   if (fabs(Point[0])>75.*mm || fabs(Point[1])>405.*mm || fabs(Point[2])>705.*mm) break;
   if (fabs(Point[0]-x[i]*mm)<5.*mm && fabs(Point[1]-y[i]*mm)<5.*mm && fabs(Point[2]-z[i]*mm)<5.*mm)
{   

    Bfield[0]=Bx[i]*gauss;
    Bfield[1]=By[i]*gauss;
    Bfield[2]=Bz[i]*gauss;
    break;}
   else if (-fabs(Point[0]-x[i]*mm)<5.*mm && fabs(Point[1]-y[i]*mm)<5.*mm && fabs(Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=-By[i]*gauss;
    Bfield[2]=-Bz[i]*gauss;
    break;}
    else if (fabs(Point[0]-x[i]*mm)<5.*mm && fabs(-Point[1]-y[i]*mm)<5.*mm && fabs(Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=-By[i]*gauss;
    Bfield[2]=Bz[i]*gauss;
    break;}
    else if (fabs(-Point[0]-x[i]*mm)<5.*mm && fabs(-Point[1]-y[i]*mm)<5.*mm && fabs(Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=By[i]*gauss;
     Bfield[2]=-Bz[i]*gauss;
    break;}
      else if (fabs(Point[0]-x[i]*mm)<5.*mm && fabs(Point[1]-y[i]*mm)<5.*mm && fabs(-Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=By[i]*gauss;
     Bfield[2]=-Bz[i]*gauss;
    break;}

     else if (fabs(-Point[0]-x[i]*mm)<5.*mm && fabs(Point[1]-y[i]*mm)<5.*mm && fabs(-Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=-By[i]*gauss;
     Bfield[2]=Bz[i]*gauss;
    break;}
    else if (fabs(Point[0]-x[i]*mm)<5.*mm && fabs(-Point[1]-y[i]*mm)<5.*mm && fabs(-Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=-By[i]*gauss;
     Bfield[2]=-Bz[i]*gauss;
    break;}
    else if (fabs(-Point[0]-x[i]*mm)<5.*mm && fabs(-Point[1]-y[i]*mm)<5.*mm && fabs(-Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=By[i]*gauss;
     Bfield[2]=Bz[i]*gauss;
    break;}
 }
}
-------------------------------------------------------------------------
   Now, I have arranged the code to open and read the files once during initialization and keep the results around to use in GetFieldValue(). The program runs much faster.
   The funny thing is: IT GOT STUCK AGAIN!
   Just like the situation before I made this improvement. I'd like to mention that the main program is fine. Because if I change the magnetic field in a simple way, the first event ends quickly and the second comes up.
   But now, there never has been an end for my first event!
   I really want to know what the matter is.

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

1 Question: Is This the Reason?   (Zukai Wang - 20 Feb, 2009)
 Add Message Add Message
to: "Re: The Strangest Problem I've ever met"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]