Geant4 Version: 11.4.0
Operating System: Ubuntu 24.04.3 LTS
Compiler/Version: gcc 13.3.0
CMake Version: 3.28.3
Dear Geant4 users,
I’ll try to explain my problem as clearly as possible, but I’m a little lost…
I’m discovering Geant4 and more precisely Geant4-DNA. I started from the example “molcounters-basic”. To the cell with its nucleus and mitochondria, I added a vesicle with in its inside a cluster of nanoparticles (NPs). The vesicle and the cluster are both placed randomly. The problem I encounter started when I tried to move the source from the center of the nucleus into the cluster of NPs. I obtained the following error message :
========================================================================================
→ G4TaskRunManager::CreateAndStartWorkers() → Creating 2 tasks with 1 events/task…G4WT0 > ### Run 0 starts on worker thread 0.
G4WT0 > ### Run 0 starts.
G4WT0 > → Event 0 starts with initial seeds (13049039,61775110).
G4WT0 > === Event 0 ===
G4WT0 > Vesicle center: (0.0021726,-0.00311178,-0.00473642)
G4WT0 > NanoCluster center: (0.00236337,-0.00285309,-0.00497114)
G4WT0 > Source local pos : (-3.77016e-05,2.90647e-05,2.44878e-05)
G4WT0 > Source global pos: (0.00232567,-0.00282403,-0.00494665)
Erreur de segmentation (core dumped)
If I understand correctly, the molcounters-basic example uses multi-threading. I tried in single-threading :
//auto* runManager = G4RunManagerFactory::CreateRunManager();
auto* runManager = new G4RunManager;
And the problem became the following :
==================================================================
G4VisManager: Using G4TrajectoryDrawByCharge as fallback trajectory model.
See commands in /vis/modeling/trajectories/ for other options.Run 0 starts.
Run 0 starts.
→ Event 0 starts.
=== Event 0 ===
Vesicle center: (-0.00130582,-0.00266929,0.0042861)
NanoCluster center: (-0.00117777,-0.00233726,0.00410789)
Source local pos : (8.73874e-05,2.93046e-05,7.1352e-06)
Source global pos: (-0.00109038,-0.00230795,0.00411502)
DNAMolecularStepByStepModel will be used-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : G4MoleculeCounterTimeComparer
issued by : G4MoleculeCounterTimeComparer::operator()
Unknown comparison type
*** Fatal Exception *** core dump ***
G4Track (0x5c2ddb65b1d0) - track ID = -52, parent ID = -3
Particle type : H3O - creator process : H2O_DNAMolecularDecay, creator model : Undefined
Kinetic energy : 38.5195 meV - Momentum direction : (5.00746e-310,-0.8651,0.45141)
Step length : 6.55902e-297 fm - total energy deposit : 7.38716 eV
Pre-step point : (-0.00118125,-0.00228721,0.004017) - Physical volume : Vesicle (G4_WATER)
- defined by : e-_G4DNAVibExcitation - step status : 4
Post-step point : (-0.00118125,-0.00228721,0.004017) - Physical volume : Vesicle (G4_WATER)- defined by : e-_G4DNAElectronSolvation - step status : 4
*** Note: Step information might not be properly updated.-------- EEEE -------- G4Exception-END --------- EEEE -------
*** G4Exception: Aborting execution ***
Abandon (core dumped)
The crash seems to be caused by the molecule counters. To confirm this I tried to comment out of the ActionInitialization.cc code each
BuildMoleculeCounters
or
BuildMultipleAndCustomMoleculeCounters
and in fact this way the problem disappeared.
Has someone any idea how to help me with this ?
Below are my DetectorConstruction.cc, PrimaryGeneratorAction.cc and MoleculeCounter.cc files.
Thank you very much!
Sarah
DetectorConstruction.cc :
G4VPhysicalVolume* DetectorConstruction::Construct()
{
CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine);
CLHEP::HepRandom::setTheSeed(time(nullptr));
auto man = G4NistManager::Instance();
auto water = man->FindOrBuildMaterial("G4_WATER");
fAGuIX = Materials::CreateAGuIX();
const G4double worldXYZ = 50 * um;
auto solidWorld = new G4Box("World", 0.5 * worldXYZ, 0.5 * worldXYZ, 0.5 * worldXYZ);
auto lvWorld = new G4LogicalVolume(solidWorld, water, "World");
auto pvWorld =
new G4PVPlacement(nullptr, G4ThreeVector(), lvWorld, "World", nullptr, false, 0, true);
const auto worldVisAtt = new G4VisAttributes(G4Colour(.5, 1.0, .5));
worldVisAtt->SetVisibility(true);
worldVisAtt->SetForceWireframe(true);
lvWorld->SetVisAttributes(worldVisAtt);
auto lvCell = ConstructCell(pvWorld);
auto lvVesicle = ConstructVesicle(lvCell);
ConstructNanoCluster(lvVesicle);
return pvWorld;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* DetectorConstruction::ConstructCell(G4VPhysicalVolume* pvWorld)
{
const G4double cellRadius = 7 * um;
const G4double nucleusRadius = 4 * um;
const G4int nMitochondria = 100;
const G4double mitoA = 0.55 * micrometer;
const G4double mitoB = 0.25 * micrometer;
const G4double mitoC = 0.90 * micrometer;
auto solidCell = new G4Orb("Cell", cellRadius);
auto lvCell = new G4LogicalVolume(
solidCell, G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"), "Cell");
auto pvCell =
new G4PVPlacement(nullptr, G4ThreeVector(), "Cell", lvCell, pvWorld, false, 0, true);
auto cellVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 1.0, 0.5));
cellVisAtt->SetVisibility(false);
lvCell->SetVisAttributes(cellVisAtt);
auto solidNucleus = new G4Orb("Nucleus", nucleusRadius);
auto lvNucleus = new G4LogicalVolume(
solidNucleus, G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"), "Nucleus");
new G4PVPlacement(nullptr, G4ThreeVector(), "Nucleus", lvNucleus, pvCell, false, 0, true);
auto solidMito = new G4Ellipsoid("Mitochondria", mitoA, mitoB, mitoC);
auto lvMito = new G4LogicalVolume(
solidMito, G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"), "Mitochondria");
for (auto i = 0; i < nMitochondria; ++i) {
G4bool overlap = true;
do {
auto u = twopi * G4UniformRand();
auto v = std::acos(2 * G4UniformRand() - 1);
auto dr = G4UniformRand() * cellRadius;
auto x = dr * std::cos(u) * std::sin(v);
auto y = dr * std::sin(u) * std::sin(v);
auto z = dr * std::cos(v);
auto pos = G4ThreeVector(x, y, z);
auto phi = G4UniformRand() * 2 * pi;
auto psi = G4UniformRand() * 2 * pi;
auto rot = new G4RotationMatrix();
rot->rotateX(psi);
rot->rotateY(phi);
auto pvMito = new G4PVPlacement(rot, pos, "Mitochondria", lvMito, pvCell, false, i, false);
overlap = pvMito->CheckOverlaps(1000, 0, false);
if (overlap) {
G4PhysicalVolumeStore::DeRegister(pvMito);
}
} while (overlap);
}
return lvCell;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* DetectorConstruction::ConstructVesicle(G4LogicalVolume* lvCell)
{
const G4double vesicleRadius = 0.5 * um;
fVesicleRadius = vesicleRadius;
auto solidVesicle = new G4Orb("Vesicle", vesicleRadius);
auto lvVesicle = new G4LogicalVolume(
solidVesicle,
G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"),
"Vesicle"
);
auto vesicleVisAtt = new G4VisAttributes(G4Colour(1.0, 1.0, 0.0));
vesicleVisAtt->SetVisibility(false);
lvVesicle->SetVisAttributes(vesicleVisAtt);
G4bool overlap = true;
G4ThreeVector vesicle_pos;
while (overlap)
{
auto u = twopi * G4UniformRand();
auto v = std::acos(2*G4UniformRand() - 1);
auto r = G4UniformRand() * 7*um;
vesicle_pos = G4ThreeVector(
r * std::cos(u) * std::sin(v),
r * std::sin(u) * std::sin(v),
r * std::cos(v)
);
G4PVPlacement* pvVesicle = new G4PVPlacement(nullptr, vesicle_pos, lvVesicle, "Vesicle", lvCell, false, 0, true);
overlap = pvVesicle->CheckOverlaps(1000, 0, false);
if (overlap)
{
G4PhysicalVolumeStore::DeRegister(pvVesicle);
}
else
{
fVesiclePosition = vesicle_pos;
}
}
return lvVesicle;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void DetectorConstruction::ConstructNanoCluster(G4LogicalVolume* lvVesicle)
{
const G4double nanoClusterRadius = 100 * nm;
fNanoClusterRadius = nanoClusterRadius;
auto solidNanoCluster = new G4Orb("NanoCluster", nanoClusterRadius);
auto lvNanoCluster = new G4LogicalVolume(solidNanoCluster, fAGuIX, "NanoCluster");
auto nanoClusterVisAtt = new G4VisAttributes(G4Colour(1.0, 1.0, 0.0,0.2));
nanoClusterVisAtt->SetForceSolid(true);
lvNanoCluster->SetVisAttributes(nanoClusterVisAtt);
double u = G4UniformRand();
double v = G4UniformRand();
double theta = std::acos(2.0*v - 1.0);
double phi = 2.0 * CLHEP::pi * u;
G4double margin = 2 * nm;
double r = fVesicleRadius - nanoClusterRadius - margin;
G4ThreeVector localPos(
r * std::sin(theta) * std::cos(phi),
r * std::sin(theta) * std::sin(phi),
r * std::cos(theta)
);
fNanoClusterPV = new G4PVPlacement(
nullptr,
localPos,
lvNanoCluster,
"NanoCluster",
lvVesicle,
false,
0,
true
);
fNanoClusterLocalPosition = localPos;
}
PrimaryGeneratorAction.cc :
PrimaryGeneratorAction::PrimaryGeneratorAction() : G4VUserPrimaryGeneratorAction()
{
G4int n_particle = 1;
fParticleGun = std::make_unique<G4ParticleGun>(n_particle);
// default particle kinematic
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* particle = particleTable->FindParticle("e-");
fParticleGun->SetParticleDefinition(particle);
fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 0.));
fParticleGun->SetParticleEnergy(100 * keV);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
// ============================================================
// 1) Récupérer la géométrie (position + rayon de la vésicule)
// ============================================================
auto det = static_cast<const DetectorConstruction*>(
G4RunManager::GetRunManager()->GetUserDetectorConstruction());
G4ThreeVector vesicleCenter = det->GetVesiclePosition(); // centre de la vésicule
G4ThreeVector nanoClusterLocal = det->GetNanoClusterLocalPosition();
G4ThreeVector clusterGlobal = vesicleCenter + nanoClusterLocal; // position globale du cluster
G4double R = det->GetNanoClusterRadius();
// ============================================================
// 2) Tirage uniforme dans le volume d’une sphère
// (distribution volumique ∝ r²)
// ============================================================
G4double u = G4UniformRand();
G4double v = G4UniformRand();
G4double w = G4UniformRand();
// pour une génération de source uniforme dans la vésicule
G4double r = R * std::cbrt(u);
G4double theta = std::acos(2.0*v - 1.0);
G4double phi = CLHEP::twopi * w;
G4ThreeVector localPos(
r * std::sin(theta) * std::cos(phi),
r * std::sin(theta) * std::sin(phi),
r * std::cos(theta)
);
// Position globale dans le repère du monde
G4ThreeVector globalPos = clusterGlobal + localPos;
G4cout << "=== Event " << anEvent->GetEventID() << " ===" << G4endl;
G4cout << "Vesicle center: " << vesicleCenter << G4endl;
G4cout << "NanoCluster center: " << clusterGlobal << G4endl;
G4cout << "Source local pos : " << localPos << G4endl;
G4cout << "Source global pos: " << globalPos << G4endl;
// ============================================================
// 2b) Écrire la position de la source dans un CSV
// ============================================================
static bool firstCall = true;
std::ofstream outFile("SourcePositions.csv", std::ios::app);
if (outFile.is_open()) {
if (firstCall) {
outFile << "EventID,X_um,Y_um,Z_um\n"; // header
firstCall = false;
}
outFile << anEvent->GetEventID() << ","
<< globalPos.x()/um << ","
<< globalPos.y()/um << ","
<< globalPos.z()/um << "\n";
outFile.close();
}
// ============================================================
// 3) Direction isotrope
// ============================================================
G4double phiDir = CLHEP::twopi * G4UniformRand();
G4double cosT = 2.0 * G4UniformRand() - 1.0;
G4double sinT = std::sqrt(1.0 - cosT*cosT);
G4ThreeVector dir(
sinT * std::cos(phiDir),
sinT * std::sin(phiDir),
cosT
);
// ============================================================
// 4) Appliquer à la source
// ============================================================
fParticleGun->SetParticlePosition(globalPos);
fParticleGun->SetParticleMomentumDirection(dir);
// ============================================================
// 5) Générer la particule
// ============================================================
fParticleGun->GeneratePrimaryVertex(anEvent);
}
And MoleculeCounter.cc :
using G4VMoleculeCounterIndex = G4VMoleculeCounter::G4VMoleculeCounterIndex;
G4ThreadLocal std::unique_ptr<G4Navigator> MoleculeCounter::fNavigator = nullptr;
MoleculeCounter::MoleculeCounter(G4String name)
: G4VUserMoleculeCounter(name, G4VMoleculeCounter::MoleculeCounterType::Other),
fIgnoreMoleculePosition(false),
fIgnoreCopyNumbers(true)
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MoleculeCounter::InitializeUser() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
std::unique_ptr<G4VMoleculeCounterIndex> MoleculeCounter::BuildIndex(const G4Track* aTrack) const
{
if (fVerbose > 1) {
G4cout << "MoleculeCounter::BuildIndex(" << aTrack->GetTrackID() << " : "
<< GetMolecule(aTrack)->GetName() << ")" << G4endl;
}
if (fIgnoreMoleculePosition) {
return std::make_unique<MoleculeCounterIndex>(
GetMolecule(aTrack)->GetMolecularConfiguration(), nullptr);
}
else {
const G4VTouchable* touchable = aTrack->GetNextTouchable();
G4TouchableHistory* touchableHistory = nullptr;
if (touchable == nullptr) touchable = aTrack->GetTouchable();
if (touchable == nullptr) {
auto touchableHandle = G4MoleculeLocator::Instance()->LocateMoleculeTrack(aTrack);
touchable = touchableHandle();
}
if (touchable == nullptr) {
G4ExceptionDescription errMsg;
errMsg << "Molecule scorer requires a valid volume pointer."
<< " G4Track->GetVolume: " << aTrack->GetVolume()
<< " G4Navigator->LocateGlobalPointAndUpdateTouchable: nullptr" << G4endl;
G4Exception("MoleculeCounter::BuildIndex", "VOL_NOT_FOUND", FatalException, errMsg);
}
// Création de la clé spatiale
// Chaque ROS est compté selon son type, le volume où il se trouve
auto index = std::make_unique<MoleculeCounterIndex>(
GetMolecule(aTrack)->GetMolecularConfiguration(), touchable);
delete touchableHistory;
return index;
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
std::unique_ptr<G4VMoleculeCounterIndex>
MoleculeCounter::BuildIndex(const G4Track* aTrack, const G4StepPoint* aStepPoint) const
{
return std::make_unique<MoleculeCounterIndex>(
GetMolecule(aTrack)->GetMolecularConfiguration(), aStepPoint->GetTouchable());
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
std::unique_ptr<G4VMoleculeCounterIndex>
MoleculeCounter::BuildSimpleIndex(const G4MolecularConfiguration* configuration) const
{
return std::make_unique<MoleculeCounterIndex>(
configuration, nullptr);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool MoleculeCounter::GetIgnoreMoleculePosition() const
{
return fIgnoreMoleculePosition;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MoleculeCounter::SetIgnoreMoleculePosition(G4bool flag)
{
fIgnoreMoleculePosition = flag;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MoleculeCounter::SetNegativeCountsAreFatal(G4bool flag)
{
G4VMoleculeCounter::SetNegativeCountsAreFatal(flag);
// this is protected in G4VMoleculeCounter and thus, must be exposed by the derived class
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......