Geant4 Version: 11.2.2
Operating System: Ubuntu
Compiler/Version: 11.4.0
CMake Version: 3.22.1
I implemented a custom discrete process in Geant4 (WimpWeakInteraction). When I run a simulation, I notice that PostStepDoIt is called twice in a row for the same step:
-
The first call has a non-zero
StepLength -
The second call has
StepLength = 0
Both calls change the particle’s direction and energy, which is not what I want.
I want the process to only modify the particle once per physical step.
I added a guard in PostStepDoIt:
```
if (steplength == 0){
aParticleChange.Initialize(track);
return &aParticleChange;}
```
This prevents the second call from affecting the particle, but I would still like to understand why Geant4 calls the process twice in the first place.
-
I am using a custom particle.
-
GetMeanFreePath()returns a finite value.
My question:
-
Why does Geant4 call
PostStepDoIttwice? -
What is the proper way to make my process execute only once per step?
```
G4VParticleChange* WimpWeakInteraction::PostStepDoIt(const G4Track& track, const G4Step& step) {
//Steplength > 0!
G4double steplength = track.GetStepLength();
/*if (steplength == 0){
aParticleChange.Initialize(track);
return &aParticleChange;}*/
// Definiere, was nach der Wechselwirkung passiert
aParticleChange.Initialize(track);
// Isotrope Streuung
const G4Material* material = track.GetMaterial();
const G4Element* element = Element(material);
G4double m_particle = MyCustomParticle::Definition()->GetPDGMass();
G4double Z = element->GetZ();
G4double A = element->GetA() / (g/mole);
G4double m_element = G4NucleiProperties::GetNuclearMass(A,Z);
G4double B = m_element / m_particle;
G4double alpha = ((B-1)/(B+1)) * ((B-1)/(B+1));
G4ThreeVector olddirection = track.GetMomentumDirection();
G4ThreeVector rand_vec = G4RandomDirection();
G4double cosc = 2 * G4UniformRand() - 1;
G4double theta_COM = std::acos(cosc);
G4double sinc = std::sin(theta_COM);
G4double theta_LAB = std::atan(sinc/(cosc+1/B));
G4double cosl = std::cos(theta_LAB);
G4double sinl = std::sin(theta_LAB);
G4ThreeVector rotaxis = rand_vec.cross(olddirection).unit();
G4ThreeVector newdirection = olddirection * cosl + rotaxis.cross(olddirection).unit() * sinl + rotaxis * (rotaxis * olddirection)*(1-cosl);
aParticleChange.ProposeMomentumDirection(newdirection);
// Ändere die Energie des WIMPs
G4double energy = track.GetKineticEnergy();
G4double newenergy = energy * (((1+alpha)+(1-alpha)*cosc)/2);
aParticleChange.ProposeEnergy(newenergy);
return &aParticleChange;
}
#include “physics.hh”
#include “particle.hh”
#include “interaction.hh”
MyPhysicsList::MyPhysicsList()
{
RegisterPhysics (new G4EmStandardPhysics());
RegisterPhysics (new G4OpticalPhysics());
RegisterPhysics (new G4DecayPhysics());
RegisterPhysics (new G4RadioactiveDecayPhysics());
}
void MyPhysicsList::ConstructParticle()
{
MyCustomParticle::Definition();
}
void MyPhysicsList::ConstructProcess()
{
AddTransportation();
G4ParticleDefinition* particle = MyCustomParticle::Definition();
G4ProcessManager* pmanager = particle->GetProcessManager();
auto weakInteraction = new WimpWeakInteraction();
pmanager->AddDiscreteProcess(weakInteraction);
}
MyPhysicsList::~MyPhysicsList()
{}
```