Documentation on Multithreading

I am raising this issue in part to encourage an update to the Application Manual and partly to document this issue for others.

Similar to the extended/medical/dna/clustering example, I have instantiated a clustering algorithm in the EventAction that I wanted to fill with ionization event-data in the SteppingAction.

Problem: This caused segfaults relating to multiple memory accesses/double-freeing of data and by printing out the memory addresses of the clustering algorithm I could confirm that the class is not instantiated on each thread as desired even though Geant4 is parallelized on Event-level.

Solution: Assuming this to be best practice, I declared the clustering algorithm instance as a private EventAction-member and wrote dedicated setter-method to register ionization sites. I eventually found out that this causes Geant4 to share one clustering instance amongst threads, which is very unhandy for this application (and not memory-save).

I find this behavior a bit unintuitive as the the setter- and getter-paradigm is very common in C++ and did not find anything about this in the documentation.

Could you post the code for the Event/Stepping Actions please, along with any code that connects them up? This would help to see what we’re missing in the docs.

Sure, thanks for reaching out. The EventAction code is something like this

#ifndef EventAction_hh
#define EventAction_hh

#include "DataTypes.hh"
#include "G4Event.hh"
#include "G4UserEventAction.hh"
#include <G4ThreeVector.hh>
#include <G4Types.hh>
#include <vector>

class MyRunAction;
class AVolClustering;

class MyEventAction : public G4UserEventAction
{
  public:
    MyEventAction(MyRunAction *);
    ~MyEventAction();

    virtual void BeginOfEventAction(const G4Event *) override;
    virtual void EndOfEventAction(const G4Event *) override;

    void RegisterIonization(G4ThreeVector location, G4double edep);

  private:
    MyRunAction *fRunAction;
    AVolClustering *fClusteringAlgo;
};

#endif

and

#include "EventAction.hh"

#include "AVolClustering.hh"
#include "CustomAnalysisManager.hh"
#include "RunAction.hh"
#include <CLHEP/Units/SystemOfUnits.h>
#include <G4ThreeVector.hh>

MyEventAction::MyEventAction(MyRunAction *runAction)
{
    fRunAction = runAction;

    double delta = 1.500888363;    // nm
    double e_min = .0;    // all energy transfers from ionizations are registered

    fClusteringAlgo = new AVolClustering(e_min, delta);
}

MyEventAction::~MyEventAction()
{
    delete fClusteringAlgo;
}

void MyEventAction::BeginOfEventAction(const G4Event *)
{
    fClusteringAlgo->Purge();
}

void MyEventAction::EndOfEventAction(const G4Event *)
{

    fClusteringAlgo->RunClustering();

    // ...

    fClusteringAlgo->Purge();
}

void MyEventAction::RegisterIonization(G4ThreeVector location, G4double edep)
{
    fClusteringAlgo->RegisterDamage(location, edep);
}

AVolClustering is a clustering algorithm that has an interface like Geant4’s DBSCAN-implementation. If it is a private member, like here, that is accessed via MyEventAction::RegisterIonization(…), Geant4 seems to share the instance fClusteringAlgorithm among threads. This causes two problems:

  1. Different workers write/delete the same memory
  2. Clustering requires ionization sites to be gathered over an event, so only ionization sites from one primary particle can form a cluster

If I daclare fClusteringAlgorithm as public instead (and remove the then unnecessary ::RegisterIonization-method), fClusteringAlgorithm is instantiated once per thread (which I confirmed by seeing different memory addresses for the instances).

The SteppingAction is simply:

#ifndef SteppingAction_hh
#define SteppingAction_hh

#include "EventAction.hh"
#include "G4Step.hh"
#include "G4UserSteppingAction.hh"
#include "G4VProcess.hh"

class G4GenericMessenger;
class MyEventAction;

class MySteppingAction : public G4UserSteppingAction
{
  public:
    MySteppingAction(MyEventAction *);
    ~MySteppingAction();

    virtual void UserSteppingAction(const G4Step *);
    // virtual void Initialize();

  private:
    MyEventAction *fEventAction;
};

#endif

and

#include "SteppingAction.hh"
#include "G4RunManager.hh"
#include <G4ThreeVector.hh>
#include <cmath>
#include <cstring>

MySteppingAction::MySteppingAction(MyEventAction *eventAction)
{
    fEventAction = eventAction;
}

MySteppingAction::~MySteppingAction()
{
}

void MySteppingAction::UserSteppingAction(const G4Step *step)
{

    G4int eventID  = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
    G4double edep  = step->GetTotalEnergyDeposit();
    G4Track *track = step->GetTrack();
    G4StepPoint *postStep = step->GetPostStepPoint();

    if ( track && track->GetNextVolume() ) {
        if ( std::strcmp(track->GetNextVolume()->GetName(), "roi_PV") == 0 ) {

            G4int processFlag =
                postStep->GetProcessDefinedStep()->GetProcessSubType();

            if ( processFlag == 12 || processFlag == 13 || processFlag == 53 ) {
                // 12 = PhotoElectric
                // 13 = Compton
                // 53 = DNAIonisation

                fEventAction->RegisterIonization(postStep->GetPosition(), edep);
            }
        }
    }
}

I can share further details/code with any mods/devs, if necessary.

Thanks! Could you also post the concrete G4VUserActionInitialization where these classes are built for the main and worker threads please?

Sure. Actually I have noted the error reappearing even though I did not set fClusteringAlgo as private. I am sorry for the confusion. In case you spot the culprit, I’d be happy to know, otherwise my plan B is to move the the clustring instance into the AnalysisManager that I wrote (inheriting from the same base classes as G4GenericAnalysisManager etc.). The G4ThreadLocalSingleton-instances have always been instantiated thread-locally.

Here’s the ActionInitialization I use:

#ifndef ActionInitialization_hh
#define ActionInitialization_hh

#include "G4VUserActionInitialization.hh"

class MyPrimaryGenerator;
class MyRunAction;
class MyEventAction;
class MySteppingAction;

class MyActionInitialization : public G4VUserActionInitialization
{
  public:
    MyActionInitialization();
    ~MyActionInitialization();

    virtual void BuildForMaster() const;
    virtual void Build() const;

    MyRunAction* GetRunAction() const;

  private:
    MyPrimaryGenerator* fGenerator;
    MyRunAction* fRunAction;
    MyEventAction* fEventAction;
    MySteppingAction* fSteppingAction;
};

#endif

and

#include "ActionInitialization.hh"
#include "EventAction.hh"
#include "Generator.hh"
#include "RunAction.hh"
#include "SteppingAction.hh"

MyActionInitialization::MyActionInitialization()
{
    fGenerator      = new MyPrimaryGenerator();
    fRunAction      = new MyRunAction();
    fEventAction    = new MyEventAction(fRunAction);
    fSteppingAction = new MySteppingAction(fEventAction);
}

MyActionInitialization::~MyActionInitialization()
{
    delete fGenerator;
    delete fRunAction;
    delete fEventAction;
    delete fSteppingAction;
}

void MyActionInitialization::BuildForMaster() const
{
    SetUserAction(fRunAction);
}

void MyActionInitialization::Build() const
{
    SetUserAction(fGenerator);
    SetUserAction(fRunAction);
    SetUserAction(fEventAction);
    SetUserAction(fSteppingAction);
}

MyRunAction* MyActionInitialization::GetRunAction() const
{
    return fRunAction;
}

Thanks for the additional code which clarifies the issue (and where we could add docs). As each action is constructed in the MyActionInitialization constructor, they are only constructed on the main thread. However, MyActionInitialization::Build() runs on each worker thread individually, and thus the actions are shared across all threads leading to the errors noted. (Note that you theoretically can do things this way, but you’d have to use mutexes in each actions member functions to avoid date races, and thus lose throughput).

Instead, to have individual actions per thread, one needs to do something like:

void MyActionInitialization::BuildForMaster() const
{
    SetUserAction(new MyRunAction);
}

void MyActionInitialization::Build() const
{
    SetUserAction(new MyPrimaryGeneratorAction);
    auto runAction = new MyRunAction;
    auto eventAction = new MyEventAction(runAction);
    auto steppingAction = new MySteppingAction(eventAction);

    SetUserAction(runAction);
    SetUserAction(eventAction);
    SetUserAction(steppingAction);
}

It’s not clear what you access the run action for, but there are methods to, e.g. merge results etc from across the worker threads.

Oh wow, that fixes it!

The MyActionInitialization::GetRunAction-method was so that I can get access RunAction-methods from within the main function which served my purpose. I have at this point already written my own AnalysisManager-class to bypass that a bit more elegantly. But it was a bit reckless to change the instantiation of the UserActions like this.

Thank you very much!

1 Like