Problems loading a simple external electric field

Geant4 Version: geant4-11-01-patch-03 [MT] (10-November-2023)
Operating System: WSL:Ubuntu

Dear forum members,

for a bit more than a month I have been trying to load an external electric field into GEANT4. I have been reading other threads on the forum on how to do it but I have come to a point where I’m stuck and have little ideas how to proceed further. Here are some of the issues I have encountered:

  1. When I try to use the command :
    /vis/scene/add/electricField i always get the “No field in this extent” message
  2. I included debugging steps and the field is being loaded correctly and the interpolation works just fine as i can print field values at certain places, and they seem to be correct.
  3. I have a source of electrons and whenever I simulate 1 particle the program either gets stuck, or the electron will accelerate the next step to something unphysical like 6^97 MeV and just spin in place. The electron has 50 keV starting velocity.

I cleaned my code and uploaded here all the files, which can be compiled to further test. It’s a simple simulation world with an electron source. Any help/tips would be appreciated!

PS: If I declare the field in G4 using the G4UniformElectricField it works without a problem.
uniform_field.txt (351.8 KB)

action.cc (268 Bytes)
action.hh (369 Bytes)
CMakeLists.txt (374 Bytes)
construction.cc (2.6 KB)
construction.hh (1.1 KB)

generator.cc (820 Bytes)
generator.hh (650 Bytes)
MyElectricField.cc (2.3 KB)
MyElectricField.hh (757 Bytes)
sim.cc (2.3 KB)

It didn’t let me upload everything so here are the physicslist declarations for completeness:
physics.cc (1.2 KB)
physics.hh (2.6 KB)

Hi Alex

Indeed, there have been some recent posts. What I would like to say, from the point of view of the visualisation, the command /vis/scene/add/electricField merely interrogates the field, and displays it. It checks that it is in the scene’s extent, i.e., within the detector - or within a user-defined extent.

If I recall, there was an issue with definition of the field. An electric field uses locations 3-5 of the field vector (roughly - I am not the expert here). The locations 0-2 are reserved for magnetic field.

John

Hi John,

thanks for your reply. I understand how this could cause visualisation problems. ,The other issue with the particle getting accelerated to unphysical values or that the program just keeps getting stuck is the big problem, and I thought these two problems might be correlated. I have a very simple implementation of the field in my constructor.cc file :

void MyDetectorConstruction::ConstructField(const G4String& fieldFilename)
{
    fElectricField = new MyElectricField(fieldFilename);
    
    G4EqMagElectricField* equation = new G4EqMagElectricField(fElectricField);
    G4MagIntegratorStepper* stepper = new G4ClassicalRK4(equation, 8);

    G4double minStep = 0.01*mm;
    G4MagInt_Driver* driver = new G4MagInt_Driver(minStep, stepper, stepper->GetNumberOfVariables());
    fChordFinder = new G4ChordFinder(driver);

    fFieldManager = new G4FieldManager(fElectricField);
    fFieldManager->SetChordFinder(fChordFinder);

    // assign the field manager to the logical volume
    logicWorld->SetFieldManager(fFieldManager, true);

    G4TransportationManager::GetTransportationManager()->GetFieldManager()->SetDetectorField(fElectricField);
    G4TransportationManager::GetTransportationManager()->GetFieldManager()->SetChordFinder(fChordFinder);
  
    // verbosity commands
    G4UImanager* UImanager = G4UImanager::GetUIpointer();
    UImanager->ApplyCommand("/tracking/verbose 1");
    UImanager->ApplyCommand("/step/verbose 1");
}  

I was wondering what could possibly go wrong here that will result in this:

* G4Track Information:   Particle = e-,   Track ID = 1,   Parent ID = 0
*********************************************************************************************************

Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
    0        0        0        0      0.05        0        0         0   physWorld initStep
    1      -60   0.0345    0.396  1.36e+98 7.26e-24     62.1      62.1  OutOfWorld Transportation

Do you have an idea? the field values are very small (1000 V/m only Ex components). I have also tried different integration steppers, and different minStep values from 1.0 mm to 0.0001 mm. Unfortunately every parameter I have tried so far didn’t seem to fix the issue.

Alex

Hi Alex

The prime suspect is MyElectricField. Have you thoroughly debugged that?

John

Hi John,

I am printing field values at random points to make sure the interpolation is fine. It does print out what the values I expect, 1000 V/m for Ex and 0 for the rest. I am doing a simple nearest-neighbor interpolation in this case. Here’s my definition of MyElectricField:

#include "MyElectricField.hh"
#include <fstream>
#include <sstream>
#include <cmath>
#include <iostream>

MyElectricField::MyElectricField(const G4String& filename)
{
    LoadFieldData(filename);
}

MyElectricField::~MyElectricField() {}

void MyElectricField::GetFieldValue(const G4double point[4], G4double* field) const
{
    InterpolateField(point, field);
}

void MyElectricField::LoadFieldData(const G4String& filename)
{
    std::ifstream infile(filename);
    if (!infile.is_open())
    {
        G4cerr << "Error opening file: " << filename << G4endl;
        return;
    }
    G4cout << "Successfully opened field data file: " << filename << G4endl;

    G4double x, y, z, Ex, Ey, Ez;
    std::string line;
    while (std::getline(infile, line))
    {
        if (line[0] == '%') continue; // skip header lines
        std::istringstream iss(line);
        if (!(iss >> x >> y >> z >> Ex >> Ey >> Ez)) { break; } // parses the line into variables, if fail, loop breaks
        fieldData.push_back({x, y, z, Ex, Ey, Ez}); // adds parsed values to the fieldData vector
    }
    infile.close();
}

void MyElectricField::InterpolateField(const G4double point[4], G4double* field) const
{
    // simple nearest-neighbor interpolation 
    G4double minDist = std::numeric_limits<G4double>::max(); //minDist is set to a large value
    const FieldPoint* nearestPoint = nullptr; // initialize pointer

    for (const auto& fp : fieldData) // loop over all field data points
    {
        // calculate the distance between the query point and the current field data point
        G4double dx = point[0] - fp.x;
        G4double dy = point[1] - fp.y;
        G4double dz = point[2] - fp.z;
        G4double dist = std::sqrt(dx*dx + dy*dy + dz*dz); // compute the euclidean distance to the input point

        if (dist < minDist)
        {
            minDist = dist;
            nearestPoint = &fp;
        }
    }

    if (nearestPoint) //if nearest point is found, the field array is set to values of the nearestPoint
    {
        field[0] = nearestPoint->Ex;
        field[1] = nearestPoint->Ey;
        field[2] = nearestPoint->Ez;
    }
    else
    {
         G4cerr << "Warning: no nearest point found." << G4endl; // if no nearestPoint found, values remains unchanged
    }
}

// debugging function to print field values at specified points
void MyElectricField::PrintFieldValues(const std::vector<G4ThreeVector>& points) const 
{
    for (const auto& point : points) 
    {
        G4double g4point[4] = { point.x(), point.y(), point.z(), 0.0 };
        G4double field[3] = { 0.0, 0.0, 0.0 };
        GetFieldValue(g4point, field);
        G4cout << "Field at (" << point.x()  << " mm, " << point.y()<< " mm, " << point.z() << " mm) : ("
               << field[0]  << " V/m, " << field[1] << " V/m, " 
               << field[2]  << " V/m)" << G4endl;
    }
}

and for completeness my corresponding .hh file:

#ifndef MYELECTRICFIELD_HH
#define MYELECTRICFIELD_HH

#include "G4ElectricField.hh"
#include "globals.hh"
#include <vector>

class MyElectricField : public G4ElectricField
{
public:
    MyElectricField(const G4String& filename);
    virtual ~MyElectricField();

    virtual void GetFieldValue(const G4double point[4], G4double* field) const override;
    // method to print field values at specified points
    void PrintFieldValues(const std::vector<G4ThreeVector>& points) const;

private:
    struct FieldPoint {
        G4double x, y, z;
        G4double Ex, Ey, Ez;
    };

    std::vector<FieldPoint> fieldData;

    void LoadFieldData(const G4String& filename);
    void InterpolateField(const G4double point[4], G4double* field) const;

};

#endif

Hi Alex

Ah! This should read

        field[0] = 0.;
        field[1] = 0.;
        field[2] = 0.;
        field[3] = nearestPoint->Ex;
        field[4] = nearestPoint->Ey;
        field[5] = nearestPoint->Ez;

with corresponding adjustment in your PrintFieldValues. Not sure where that is documented, but it’s clear from the vis modelling function

void G4ElectricFieldModel::
GetFieldAtLocation(const G4Field* field, const G4Point3D& position,
                   G4double time, G4Point3D& result) const {
  if (!field) return;                   // No action if no field

  G4double xyzt[4] = { position.x(), position.y(), position.z(), time };
  G4double BEvals[6] = {0.};            // Field returns {Bx,By,Bz,Ex,Ey,Ez}
  field->GetFieldValue(xyzt, BEvals);

  result.set(BEvals[3], BEvals[4], BEvals[5]);
  return;
}

Give it a try!

John

1 Like

Wow! I can’t thank you enough!!
It seems that the program is working now. The visualisation command works fine as well and it doesn’t throw me the error of “No field in this extent”.

Can you please explain me what was wrong in the code? and how this fixed it?

        field[0] = 0.;
        field[1] = 0.;
        field[2] = 0.;
        field[3] = nearestPoint->Ex;
        field[4] = nearestPoint->Ey;
        field[5] = nearestPoint->Ez;

I’ve been only using G4 for a few months now and any explanation is very much appreciated. Also if I would want to load a magnetic field, would I have to use the same approach?

Hi Alex

That’s great. The first 3 components of field are the magnetic field.

I’m sure this is all documented somewhere, or illustrated in an example, but I leave that to you to discover. I just looked at how the vis modelled it.

Good luck

John