Creating non-uniform electric field

Dear John,
I have tried running the program with a uniform electric field and I still seem to have the same problem. I have uploaded an image of the call stack - if its of any use (I can´t seem to copy and paste the text). I note that the problem does not seem to occur after I type
/run/beamOn 1
for the first time. It only happens if I type
/run/beamOn 1
a second time

Best Wishes
Chris

Hi Chris

Looks like a problem at run initialisation, when the run manager attempts to delete an event. Why should the run manager be deleting an event at run intialisation? I think we need to alert our run and tracking experts on this one.

Please confirm it happens only for smooth trajectories.

Also look out for warning messages - maybe way back.

Also, any chance you build in debug mode so we get a more complete traceback:
-DCMAKE_BUILD_TYPE=Debug

John

Dear John,
I can confirm that it seems to happen only for smooth trajectories.

I have the following warning message right at the start. Sorry, I should have noticed it before. I will look into how to add an exception handler.

-------- WWWW ------- G4Exception-START -------- WWWW -------

*** ExceptionHandler is not defined ***
*** G4Exception : UI0002
issued by : G4UIExecutive::G4UIExecutive()
Specified session type is not build in your system,
or no session type is specified.
A fallback session type is used.
*** This is just a warning message. ***
-------- WWWW ------- G4Exception-END -------- WWWW -------

I have recompiled in Debug and I attach the screenshots of various points in the call stack.
Best Wishes
Chris

Dear John,
Looking into the problem a little more, I see that the program only crashes after an electron has interacted with the sample. If I fire an electron away from the sample, then I can repeatedly do this without the program crashing. This makes me think that something in the physics code (which I have inherited) is causing the problem. Most likely some area of memory is being overwritten which happens to be something to do with the the G4SmoothTrajectoryPoint. Perhaps setting a watch on certain memory locations might help. Would it be possible to use valgrind ? I could try tonight on my Linux system…

i haven´t figured out how to stop the warning about G4UIExecutive and the Exception handler not defined, but it does day a fallback session type is used…

Best Wishes
Chris

Hi Chris

Have you made any progress with this?

The G4Exception is, as it says, just a warning message. What session type, if any, are you specifying? Early on there should be a line that says something like
Available UI session types: [ Qt, Xm, GAG, tcsh, csh ]
(See Application Developers Guide, How to Select Interface in Your Applications)

You may specify a session type with environment variables, e.g.,
export G4UI_USE_QT=1
or for csh
setenv G4UI_USE_QT 1
and if not available, it will choose the next most highly featured.

You may also specify the session type with ~/.g4session, a file in your home directory that contains something like this (this is the first few lines of mine):

Qt # Default session
exampleB1 tcsh
#exampleB1 GAG
exampleB2a tcsh
exampleB2b Xm

Hope that helps.

John

Hi John,

I am on holiday at the moment, but I will try what you suggest when I get back. I didn't make so much progress with it. I discovered something in my program that says something along the lines of "conditional jump depends on uninitialised variable" and gives a line and a file. The file has nothing to do with visualisation and although there was an if statement at the given line, when I tried to print out the variables they either did not print out, or they were sensible numbers. I am not certain that this is the cause of the crash when trajectory drawing. I decided to concentrate on something else for the moment as I only need the visualisation to check that everything is working correctly and not for doing the actual simulations and it only crashes after I do a run/beamOn command for the second time.- although it is rather disconcerting to have a bug in the program that causes a crash. 

I have hit a slight snag with the non-linear electric field, but we can get around it by launching electrons  within COMSOL and then launching the electrons in Geant4 once they hit the sample surface.

Best Wishes

Chris

Dear John,
Sorry I didn´t get back to you sooner, but I have had some experimental and design work to do in addition to the Geant4 stuff. The line about UI session types says:

Available UI session types: [ Win32, GAG, csh ]

on the windows PC. I can check what it says on the Linux later.

A couple of problems I have are that

  1. I am using a 3D grid of points to represent the electric field. As the field considerably varies in strength depending on how far away the tip is, I use different grids at different scales. However, when I launch the electrons from the tip, the surface of the tip does not coincide with any grid points. Hence the electric field near the tip is interpolated incorrectly as the field inside the tip is zero. I can get around this for the moment by launching electrons from the sample surface. If I need to launch electrons from the tip, the current plan is to launch using COMSOL and note where the electrons strike the sample surface and then launch electrons from the positions as stored by COMSOL.
  2. As a test, I thought I would launch electrons using COMSOL and Geant4 from the same point on the surface with the same conditions. I should see the same results (more or less) for electrons flying through the vacuum. However, at the moment I do not. Secondary electrons generated within the material seem to fly in arcs if they escape into the vacuum, but primary electrons do not seem to feel the electric field and fly in straight lines right out of the tip-sample region (even for very low energies).
  3. I still have the problem with the program crashing. I just exit the program each time if I want to view a new set of electron trajectories.

Best Wishes
Chris

Did you implement your 3D grid with an interpolation function? Or are you returning the value at each coordinate, i.e., zero in between grid point?

Hi Mike,
I am implementing a 3D interpolation function which I have copied from the purging magnet example. In addition, I have implemented a messenger function /field/getField so I can determine the field at any given point.

At the bottom of this message is the code of the interpolation function (if your´re interested !). I wrote some code so that I could have a coarser field further away from the tip, but this adds a complication, so I have returned to the simpler code using just one scale for the interpolation.

I have also attached a screen shot of the trajectories. The electrons are launched just above the surface with an energy of 1 eV. Some electrons feel the field and bend back to the surface as they should, but others just fly away as if they do not feel the field at all.

Best Wishes
Chris

void PurgMagTabulatedEField3D::GetFieldValue(const double point[4],
double *BEfield ) const
{

BEfield[0] = 0.0;
BEfield[1] = 0.0;
BEfield[2] = 0.0;

double x = point[0];
double y = point[1];
double z = point[2] + fZoffset;

#ifdef DEBUG_INTERPOLATING_FIELD
G4cout << "x, y, z: " << x << " " << y << " " << z << endl;
G4cout << "minx, miny, minz: " << minx << " " << miny << " " << minz << endl;
G4cout << "maxx, maxy, maxz: " << maxx << " " << maxy << " " << maxz << endl;
#endif

// Check that the point is within the defined region
if ( x>=minx && x<=maxx &&
y>=miny && y<=maxy &&
z>=minz && z<=maxz ) {

// Position of given point within region, normalized to the range
// [0,1]
double xfraction = (x - minx) / dx;
double yfraction = (y - miny) / dy; 
double zfraction = (z - minz) / dz;

if (invertX) { xfraction = 1 - xfraction;}
if (invertY) { yfraction = 1 - yfraction;}
if (invertZ) { zfraction = 1 - zfraction;}

// Need addresses of these to pass to modf below.
// modf uses its second argument as an OUTPUT argument.
double xdindex, ydindex, zdindex;

// Position of the point within the cuboid defined by the
// nearest surrounding tabulated points
double xlocal = ( std::modf(xfraction*(nx-1), &xdindex));
double ylocal = ( std::modf(yfraction*(ny-1), &ydindex));
double zlocal = ( std::modf(zfraction*(nz-1), &zdindex));

// The indices of the nearest tabulated point whose coordinates
// are all less than those of the given point
int xindex = static_cast<int>(xdindex);
int yindex = static_cast<int>(ydindex);
int zindex = static_cast<int>(zdindex);

#ifdef DEBUG_INTERPOLATING_FIELD
G4cout << "Fraction x,y,z: " << xfraction << " " << yfraction << " " << zfraction << endl;
G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << endl;
G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << endl;
double valx0z0, mulx0z0, valx1z0, mulx1z0;
double valx0z1, mulx0z1, valx1z1, mulx1z1;
valx0z0= xField[xindex ][0][zindex]; mulx0z0= (1-xlocal) * (1-zlocal);
valx1z0= xField[xindex+1][0][zindex]; mulx1z0= xlocal * (1-zlocal);
valx0z1= xField[xindex ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
valx1z1= xField[xindex+1][0][zindex+1]; mulx1z1= xlocal * zlocal;
#endif

    // Full 3-dimensional version
BEfield[3] =
  xField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
  xField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
  xField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
  xField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
  xField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
  xField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
  xField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
  xField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
BEfield[4] =
  yField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
  yField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
  yField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
  yField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
  yField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
  yField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
  yField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
  yField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
BEfield[5] =
  zField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
  zField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
  zField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
  zField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
  zField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
  zField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
  zField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
  zField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;

#ifdef DEBUG_INTERPOLATING_FIELD
G4cout << "minx: " << minx << " maxx: " << maxx << " x: " << x << " BEfield[3]: " << BEfield[3] << G4endl;
G4cout << "miny: " << miny << " maxy: " << maxy << " y: " << y << " BEfield[4]: " << BEfield[4] << G4endl;
G4cout << "minz: " << minz << " maxz: " << maxz << " z: " << z << " BEfield[5]: " << BEfield[5] << G4endl;
#endif

} else {

BEfield[3] = 0.0;
BEfield[4] = 0.0;
BEfield[5] = 0.0;

#ifdef DEBUG_INTERPOLATING_FIELD
// G4cout << "minx: " << minx << " maxx: " << maxx << " x: " << x << G4endl;
// G4cout << "miny: " << miny << " maxy: " << maxy << " y: " << y << G4endl;
// G4cout << "minz: " << minz << " maxz: " << maxz << " z: " << z << G4endl;
#endif

}
}

In Geant4, 1 eV electrons will either get stopped instantly in material (full energy loss in one step), or will travel indefinitely in vacuum (which I’m guessing is what you want and have set up).

One possibility may be that because there are no interactions, G4Transportation is making a single step from the starting point out to the volume (world?) boundary, so it doesn’t see the field variations. You might need to add G4StepLimiterPhysics to your physics list, and specify a maximum step length that is comparable to your field variations.

Hi Mike,
Thanks for the reply. I implemented what you suggested. My PhysicsList inherited from G4VUserPhysicsList, so I had to change it so that it inherited from G4VModularPhysicsList. However, it seems there is a problem with my GetFieldValue function as it returns zero for all points on the z axis (which it shouldn´t) and it just so happened that I was launching my electrons from various points on the z axis. Things behave more sensibly if I launch at other points.
If John Allinson is reading this, then in my efforts to find out what was going on, I expanded the size of the problem to the metre scale (instead of nm) and I changed the PhysicsList to QBBC as opposed to my own home grown PhysicsList and I increased the simulation energies. I also was having a problem with the program crashing when I visualised the trajectories for a second time. This doesn´t happen with the QBBC PhysicsList, so I suppose it is something about the PhysicsList I am using which is causing problems with the visualisation.

Best Wishes
Chris

Thanks for the feedback, Chris. Glad you’re up and running sensibly now.

Mike’s idea of limiting the step length might still be worth trying.

John

Hi John and Mike,
I will give that a go. One other thing cropped up while I was looking into this. The field near (i.e. up to 10 nm away) the tip is very strong, however further away (i.e. mm distances) the field can be as much as 6 orders of magnitude smaller. I am intending to generate data at fine scales near the tip and at a much coarser scale further from the tip. However, the miss distance (and maybe other parameters) also needs to be changed depending on the strength of the field. Since the miss distance is associated with a field (either electric or magnetic), I presume I will need to setup several volumes and a different field (and miss distance) in each volume ?

Best Wishes
Chris

I pass on this question. :slight_smile:

I have not used varying miss distance with the fields; I"ve only used the default value.

OK, Thanks for replies, I am exploring whether I can do the simulation with just one volume. If it doesn´t work, I will implement several volumes for the electric field.

Previously I had mentioned that my program crashes after I visualise trajectories for a second time (i.e. I do run/beamOn and I see trajectories, then I do run/beamOn again and the program crashes). This still happens and I now realise this is something to do with my class that handles elastic scattering and the home grown version of G4Trajectory. These two classes need to speak to each other, and I note there have been some significant changes to G4Trajectory in version 10.6. Hence, I looked through the code in my home grown version of G4Trajectory and then attempted to implement the modifications in that class into the new geant4 version (10.6) of G4Trajectory and recompile. Unfortunately, the program still crashes. I think it has to do with the way the "Safety" is calculated, and this seems to be calculated in a different way in 10.6. It ultimately crashes when it tries to delete fAuxiliaryPointVector in G4SmoothTrajectoryPoint, but the cause is certainly elsewhere. At any rate, I thought I had updated my code to version 10.6, but evidently this is not completely true. Presently, I am not sure how to update the code as I don´t know enough about how the safety is calculated. I can still carry on with my program as long as I don´t view the trajectories twice with run/beamOn. However, it would be nice to fix this. Would posting the output from valgrind help or does anyone have a suggestion as to how I might pursue this ?

Thanks in advance
Chris