Custom user-defined ion with decay table never decays during tracking

norach1991 · March 12, 2026, 7:48pm

Summary

My goal is to create ions with user-defined static properties and a user-defined decay scheme.
The new ion should transport and interact with matter like any other ion defined in Geant4.

Problem:
The ion never decays during tracking, even though its lifetime is extremely short.

—————————————————————————————————————————–-
Geant4 version: 11.3.2
Description: Ubuntu 24.04.2 LTS
Linux LT-25-0031 6.6.87.2-microsoft-standard-WSL2 #1 SMP PREEMPT_DYNAMIC Thu Jun 5 18:30:46 UTC 2025 x86_64 x86_64 x86_64 GNU/Linux
g++ (Ubuntu 13.3.0-6ubuntu2~24.04.1) 13.3.0
cmake version 3.28.3
——————————————————————————————————————————–

Greetings to the community,
I would like to define ions with the following parameters completely defined by the user.
G4Ions(
const G4String& aName, G4double mass,
G4double width, G4double charge,
G4int iSpin, G4int iParity,
G4int iConjugation, G4int iIsospin,
G4int iIsospin3, G4int gParity,
const G4String& pType, G4int lepton,
G4int baryon, G4int encoding,
G4bool stable, G4double lifetime,
G4DecayTable *decaytable, G4bool shortlived,
const G4String& subType =“”,
G4int anti_encoding =0,
G4double excitation = 0.0
);
including a user-defined decay scheme.
These ions should behave exactly like any other ion implemented in Geant4 (transportation, interaction with matter, etc.).

Below are the steps I followed.

1) I defined a custom ion type inheriting from G4Ions
#ifndef DRAGONION_HH
#define DRAGONION_HH

#include “G4Ions.hh”

namespace DRAGON
{
class DRAGONIon : public G4Ions {
public:
DRAGONIon(const G4String& name,
G4double mass,
G4double width,
G4double charge,
G4int iSpin,
G4int iParity,
G4int iConjugation,
G4int iIsospin,
G4int iIsospin3,
G4int gParity,
const G4String& type,
G4int lepton,
G4int baryon,
G4int encoding,
G4bool stable,
G4double lifetime,
G4DecayTable* decaytable,
G4bool shortlived,
G4int Z,
G4int A);
virtual ~DRAGONIon() = default;
};
}
#endif

#include “DRAGONIon.hh”

namespace DRAGON
{
DRAGONIon::DRAGONIon(const G4String& name,
G4double mass,
G4double width,
G4double charge,
G4int iSpin,
G4int iParity,
G4int iConjugation,
G4int iIsospin,
G4int iIsospin3,
G4int gParity,
const G4String& type,
G4int lepton,
G4int baryon,
G4int encoding,
G4bool stable,
G4double lifetime,
G4DecayTable* decaytable,
G4bool shortlived,
G4int Z,
G4int A)
: G4Ions(name, mass, width, charge, iSpin, iParity, iConjugation,
iIsospin, iIsospin3, gParity, type, lepton, baryon,
encoding, stable, lifetime, decaytable, shortlived)
{
SetAtomicNumber(Z);
SetAtomicMass(A);
}
}

2) Inside a method called ureact() (which is called from ConstructParticle()), I create the ions involved
//C
//C → Define Ne19++++ ion
//C

ipart = 61;
aamass = 19.*0.93149432;
ubuf[0] = 4;
tlif = 1000.;

if(!G4ParticleTable::GetParticleTable()->FindParticle("Ne19++++ "))
{
DRAGONIon* ion = new DRAGONIon(
"Ne19++++ ",
aamass*GeV,
0.0,
ubuf[0]eplus,
1, +1, 0, 0, 0, 0,
“nucleus”,
0, +19, ipart,
false,
tlifs,
nullptr,
false,
10,
19
);
}

3) I define and assign the decay scheme with something like this
brat[0] = 95.;
mode[0] = 84;
brat[1] = 5.;
mode[1] = 83;

G4DecayTable* decayTable2 = new G4DecayTable();
for (int i = 0; i < 10; i++)
{
if(mode[i])
{
G4VDecayChannel* decayChannel = new G4PhaseSpaceDecayChannel(
G4ParticleTable::GetParticleTable()->FindParticle(82)->GetParticleName(),
brat[i]/100.,
1,
G4ParticleTable::GetParticleTable()->FindParticle(mode[i])->GetParticleName()
);
decayTable2->Insert(decayChannel);
}
}

DRAGONIon* part1 =
dynamic_cast<DRAGONIon*>(G4ParticleTable::GetParticleTable()->FindParticle(82));
part1->SetDecayTable(decayTable2);
decayTable2->DumpInfo();

4) Immediately after that, I print the ion information and obtain
===== Ion Properties =====
Name: 16O_3
Type: nucleus
SubType:
PDG Encoding: 89
Mass: 14.9061 GeV
Charge: 6 e
Spin: 0
Magnetic Moment: 0
Parity: 1
Isospin: 0
Isospin3: 0
Stable: No
Lifetime: 4.7e-15 s
Atomic Number (Z): 8
Atomic Mass (A): 16
Excitation Energy: 0 eV

G4DecayTable: 16O_3
0: BR: 0.9997 [Phase Space] : 16O_
1: BR: 0.00027 [Phase Space] : 16O_1
2: BR: 8e-05 [Phase Space] : 16O_2

However, during tracking I would expect these ions to decay almost immediately, given their extremely short lifetime. Nevertheless, they never decay.

5) This is the physics list implementation

DRAGONPhysicsList::~DRAGONPhysicsList()
{
delete fMessenger;
delete fDecPhysicsList;
delete fEmPhysicsList;
delete fIonPhysicsList;
delete fStepLimiterPhysics;
delete fHadElastPhysics;
delete fHadPhysicsFTFP_BERT;
}

void DRAGONPhysicsList::ConstructParticle()
{
// Particle table initialization
fEmPhysicsList->ConstructParticle();
fDecPhysicsList->ConstructParticle();
fIonPhysicsList->ConstructParticle();
fHadElastPhysics->ConstructParticle();
fHadPhysicsFTFP_BERT->ConstructParticle();
// Define radioactive ion reactions
ureact();
}

void DRAGONPhysicsList::ConstructProcess()
{
G4cout << “From Physics ConstructProcess()” << G4endl;
AddTransportation();
fEmPhysicsList->ConstructProcess();
fDecPhysicsList->ConstructProcess();
fIonPhysicsList->ConstructProcess();
fStepLimiterPhysics->ConstructProcess();
fHadElastPhysics->ConstructProcess();
fHadPhysicsFTFP_BERT->ConstructProcess();
AddIonGasModels();
}

Any help or suggestions would be greatly appreciated.

Best regards.

pico · March 16, 2026, 2:45pm

Dear @norach1991 ,

What is the purpose of these DRAGONIONS ? There are ways to use a “regular” Ne-19 ion and simply hijack it’s decay table. You can do so by “overriding” the file to be used as the decay table. Whether a short [too] short lived particle is transported or decayed immediately is another question.

/process/had/rdm/setRadioactiveDecayFile Z A FILE

you can copy the file from the regular dataset, just grab z7.a19, modify it and point to it.

/Pico

norach1991 · May 11, 2026, 8:52pm

Thanks for your reply. The issue is that I need full freedom to define the values of all these parameters that define the particle and are given to the G4ParticleDefinition class’ constructor. Defining an ion via G4ParticleDefinition *na_ion = ionTable->GetIon(Z, A, 0); is not a good option because it does not offer me the level of freedom I need, so I need to define my own ions.

I call this method from ConstructParticle():

void DRAGONPhysicsList::CreateExoticIon() {
  if (!G4ParticleTable::GetParticleTable()->FindParticle("ExoticIon")) {

    G4String name = "ExoticIon";
    G4double mass = 50.0*GeV;
    G4double width = 0.0*MeV;
    G4double charge = 20.0*eplus;
    G4int spin = 0;
    G4int parity = +1;
    G4int CConjugation = 0;
    G4int isospin = 0;
    G4int isospin3 = 0;
    G4int GParity = 0;
    G4String type = "nucleus";
    G4int leptonNumber = 0;
    G4int baryonNumber = 50;
    G4int encoding = 12345;

    G4ParticleDefinition* ion = new G4ParticleDefinition(
        name, mass, width, charge, spin, parity, CConjugation,
        isospin, isospin3, GParity, type, leptonNumber, baryonNumber,
        encoding, false, 0.0, nullptr, false
    );

    G4VDecayChannel* gammaChannel1 =
        new G4PhaseSpaceDecayChannel(ion->GetParticleName(), 0.5, 2, "gamma", "16O_3");

    G4VDecayChannel* gammaChannel2 =
        new G4PhaseSpaceDecayChannel(ion->GetParticleName(), 0.5, 2, "gamma", "16O_2");

    gammaChannel1->SetVerboseLevel(3);
    gammaChannel2->SetVerboseLevel(3);

    G4DecayTable* decayTable = new G4DecayTable();
    decayTable->Insert(gammaChannel1);
    decayTable->Insert(gammaChannel2);

    ion->SetDecayTable(decayTable);
    decayTable->DumpInfo();

    G4IonTable::GetIonTable()->Insert(ion);
  }
}

And this one from ConstructProcess():

void DRAGONPhysicsList::AddExoticIonProcesses() {

  G4ParticleDefinition* ion =
      G4ParticleTable::GetParticleTable()->FindParticle("ExoticIon");

  if (!ion) {
    G4cout << "ExoticIon not found!" << G4endl;
    return;
  }

  G4ProcessManager* pm = ion->GetProcessManager();

  if (!pm) {
    pm = new G4ProcessManager(ion);
    ion->SetProcessManager(pm);
  }

  // ---- Multiple scattering ----
  G4hMultipleScattering* msc = new G4hMultipleScattering();
  pm->AddProcess(msc, -1, 1, 1);

  // ---- Ionisation ----
  G4ionIonisation* ionIoni = new G4ionIonisation();
  pm->AddProcess(ionIoni, -1, 2, 2);

  // ---- Nuclear stopping ----
  G4NuclearStopping* nucStop = new G4NuclearStopping();
  pm->AddProcess(nucStop, -1, 3, -1);

  G4cout << "Custom ion processes added to ExoticIon" << G4endl;
}

These are the only important processes I need for this user-defined ion in my simulation.
During the simulation I see:

Particle ExoticIon has a strange PDGEncoding
-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : PART102
issued by : G4ParticleDefintion::G4ParticleDefintion
Strange PDGEncoding
*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------

G4DecayTable:  ExoticIon
0:  BR:  0.5  [Phase Space]   :   gamma 16O_3
1:  BR:  0.5  [Phase Space]   :   gamma 16O_2

G4WT0 > G4PhaseSpaceDecayChannel::DecayIt()
G4WT0 > G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
G4WT0 > G4PhaseSpaceDecayChannel::TwoBodyDecayIt() - Create decay products in rest frame
G4WT0 >  ----- List of DecayProducts  -----
G4WT0 >  ------ Parent Particle ----------
G4WT0 >  Particle type - ExoticIon
G4WT0 >    mass:        50[GeV]
G4WT0 >    charge:      2[e]
G4WT0 >    Direction x: 0, y: 0, z: 0
G4WT0 >    Total Momentum = 0[GeV]
G4WT0 >    Momentum: 0[GeV], y: 0[GeV], z: 0[GeV]
G4WT0 >    Total Energy   = 50[GeV]
G4WT0 >    Kinetic Energy = 0[GeV]
G4WT0 >  MagneticMoment  [MeV/T]: 0
G4WT0 >    ProperTime     = 0[ns]
G4WT0 >  ------ Daughter Particles  ------
G4WT0 >  ----------1 -------------
G4WT0 >  Particle type - gamma
G4WT0 >    mass:        0[GeV]
G4WT0 >    charge:      0[e]
G4WT0 >    Direction x: -0.393984, y: -0.0868169, z: -0.915008
G4WT0 >    Total Momentum = 22.7781[GeV]
G4WT0 >    Momentum: -8.9742[GeV], y: -1.97752[GeV], z: -20.8421[GeV]
G4WT0 >    Total Energy   = 22.7781[GeV]
G4WT0 >    Kinetic Energy = 22.7781[GeV]
G4WT0 >  MagneticMoment  [MeV/T]: 0
G4WT0 >    ProperTime     = 0[ns]
G4WT0 >  ----------2 -------------
G4WT0 >  Particle type - 16O_3
G4WT0 >    mass:        14.9061[GeV]
G4WT0 >    charge:      6[e]
G4WT0 >    Direction x: 0.393984, y: 0.0868169, z: 0.915008
G4WT0 >    Total Momentum = 22.7781[GeV]
G4WT0 >    Momentum: 8.9742[GeV], y: 1.97752[GeV], z: 20.8421[GeV]
G4WT0 >    Total Energy   = 27.2219[GeV]
G4WT0 >    Kinetic Energy = 12.3158[GeV]
G4WT0 >  MagneticMoment  [MeV/T]: 0
G4WT0 >    ProperTime     = 0[ns]
G4WT0 >  ----- End List of DecayProducts  -----

From UserSteppingAction():

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

  const G4Track* track5 = step->GetTrack();
  G4ParticleDefinition* particle5 = track5->GetDefinition();

  if (particle5) {
    G4ProcessManager* processManager = particle5->GetProcessManager();
    if (processManager) {
      G4ProcessVector* processVector = processManager->GetProcessList();
      if (processVector) {

        G4cout << "Processes for particle '"
               << particle5->GetParticleName() << "':" << G4endl;

        for (size_t i = 0; i < processVector->size(); i++) {
          G4VProcess* process = (*processVector)[i];
          if (process) {
            G4cout << "  - Process " << i << ": "
                   << process->GetProcessName() << G4endl;
          }
        }
      }
    }
  }
}

And during simulation I get:

Processes for particle 'ExoticIon':
  - Process 0: Transportation
  - Process 1: Decay
  - Process 2: StepLimiter
  - Process 3: UserSpecialCut
  - Process 4: msc
  - Process 5: ionIoni
  - Process 6: nuclearStopping

MY CONCERN:
I was wondering if the way I proceeded is correct and works even despite the issue that the PDG encoding does not follow the Geant4 internal structure. This might be a problem for the evaluation of the A and Z values for this ion (important for nuclear processes like nuclear stopping and others).

Thank you!

pico · May 29, 2026, 3:00pm

Hi,

this is not something I have tried myself. With very low priority I can have a look at what you are attempting. If in the meantime you work it out, please let us know.

/Pico