Physics list for decay chain and gamma cascades

Getting Started
1

Please fill out the following information to help in answering your question, and also see tips for posting code snippets. If you don’t provide this information it will take more time to help with your problem!

_Geant4 Version:_11.2.2
_Operating System:_Ubuntu 24.04.1
_Compiler/Version:_13.3.0
_CMake Version:_3.28.3

Hi everyone,

I am trying to simulate gamma rays (< 5 MeV) from the decay chain of nuclides like U-238. I have detectors that detect the emitted gamma rays. Both the timing (mean life of excited states in gamma cascades) and angles of emitted gamma rays are important for my application. X-rays from decays are also important.

I have looked at examples like rdecay01 and slides from previous Geant4 courses to write my code. I am using G4EmStandardPhysics_option4.hh because I found it includes production of X-rays.

My problem: when I am trying to set deex->SetCorrelatedGamma(true);, I get a segmentation fault. If I comment that line out it works fine.

  1. Could someone with more experience first tell me if my physics lists are reasonable for my application? I have included G4HadronPhysicsFTFP_BERT.hh, but do I need it?
  2. What could be the reason for the segmentation fault? Am I missing an #include or similar?

Thank you for your time.

This is my code:
EPhysicslist.cc

#include "EPhysicsList.hh"

EPhysicsList::EPhysicsList()
{
    // register the physics lists with G4VModularPhysicsList
    RegisterPhysics(new G4EmStandardPhysics_option4());
    RegisterPhysics(new G4DecayPhysics());
    RegisterPhysics(new G4RadioactiveDecayPhysics());
    RegisterPhysics(new G4HadronPhysicsFTFP_BERT());

    // update G4NuclideTable time limit
    const G4double meanLife = 1*picosecond;
    G4NuclideTable::GetInstance()->SetMeanLifeThreshold(meanLife);
    G4NuclideTable::GetInstance()->SetLevelTolerance(1.0*eV);

    // define flags for nuclear gamma de-excitation model
    G4DeexPrecoParameters* deex = G4NuclearLevelData::GetInstance()->GetParameters();
    deex->SetCorrelatedGamma(true);
    deex->SetStoreAllLevels(true);
    deex->SetInternalConversionFlag(true);
    deex->SetIsomerProduction(true);
    deex->SetMaxLifeTime(meanLife);
}

EPhysicsList::~EPhysicsList()
{
    
}

and
EPhysicslist.hh

#ifndef EPHYSICSLIST_HH
#define EPHYSICSLIST_HH

#include "G4VModularPhysicsList.hh"

#include "G4EmStandardPhysics_option4.hh" 
#include "G4DecayPhysics.hh"
#include "G4RadioactiveDecayPhysics.hh"
#include "G4HadronPhysicsFTFP_BERT.hh"

#include "G4NuclideTable.hh"
#include "G4DeexPrecoParameters.hh"
#include "G4NuclearLevelData.hh"

#include "G4SystemOfUnits.hh"

class EPhysicsList : public G4VModularPhysicsList
{
public:
    EPhysicsList();
    ~EPhysicsList();
};

#endif
2

Your physics list looks reasonable. It’s not at all obvious from the code why you’d get a segfault at the line you claim. Run in the debugger, and when it segfaults, get a backtrace to see exactly what is failing where.

It may be that you should probably be doing all of that configuration work in ConstructProcess(), not in the constructor.

I would encourage you to follow C++ best practices, and move all of your #includes except for the base class to the .cc file. The class definition (.hh file) doesn’t need to know any of those other classes. It may not be “important” here, but having good habits will serve you much better as you develop more complex code.

1 Like
3

Something you’ll want to keep in mind about this. Geant4 does not model secular equilibrium.

If you use G4ParticleGun (or GPS, or whatever), and start with a U-238 nucleus at rest, each event will contain the complete decay chain, one nucleus each, all the way down to lead. The daughter gammas, alphas, etc. will have GlobalTime() values of a few E+25 ns (in G4 units), because each one will be assigned a random time based on the half-life of the decaying nucleus.

If you are wanting to model what a detector will see from a sample of purified uranium, then you don’t want the decay chain, you just want the U-238 decay itself, and then model whatever induced interactions happen in the material, container, etc.

If you’re modelling rock, such as the walls of an underground facility, then you’ll want to use software (such as Sources4C or RadSrc) to solve the Bateman equations for secular equilibrium (or just look up the result), and generate a gamma spectrum based on the relative populations of the daughters.

1 Like
4

Thank you for the feedback! I have spent some time trying to figure this segfault out, but I still need to ask for help I think.

Now I have discovered that most nuclides work fine with deex->SetCorrelatedGamma(true);. Some nuclides are problematic and result in a segfault, like Bi-214 for some reason. I set nuclides (one decay at a time with nucleusLimits) in a macro file:

/run/numberOfThreads 6

/run/initialize

/process/had/rdm/thresholdForVeryLongDecayTime 1.0e+60 year

/gun/particle ion
/gun/ion 83 214 0 0
/process/had/rdm/nucleusLimits 214 214 83 83

/run/printProgress 1000
/run/beamOn 10000

I have tried to use gdb and backtrace full but struggle to get from the output to the reason for the segfault. Just to share what I did:

  1. gdb ./sim (my Geant4 code main is in sim.cc)
  2. run run.mac (The macro above is called run.mac)
  3. It encounters the segfault.
  4. backtrace full
  5. Output:
#0  0x00007ffff65bb938 in G4PolarizationTransition::SampleGammaTransition(G4NuclearPolarization*, int, int, int, int, double, double&, double&) ()
   from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4processes.so
No symbol table info available.
#1  0x00007ffff65b3864 in G4GammaTransition::SampleDirection(G4Fragment*, double, int, int, int) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4processes.so
No symbol table info available.
#2  0x00007ffff65b364e in G4GammaTransition::SampleTransition(G4Fragment*, double, double, int, int, int, int, bool, bool) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4processes.so
No symbol table info available.
#3  0x00007ffff65b65bb in G4PhotonEvaporation::GenerateGamma(G4Fragment*) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4processes.so
No symbol table info available.
#4  0x00007ffff65b7ede in G4PhotonEvaporation::EmittedFragment(G4Fragment*) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4processes.so
No symbol table info available.
#5  0x00007ffff69fd0ce in G4ITDecay::DecayIt(double) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4processes.so
No symbol table info available.
#6  0x00007ffff6a0321a in G4RadioactiveDecay::DoDecay(G4ParticleDefinition const&, G4DecayTable*) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4processes.so
No symbol table info available.
#7  0x00007ffff6a038de in G4RadioactiveDecay::DecayAnalog(G4Track const&, G4DecayTable*) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4processes.so
No symbol table info available.
#8  0x00007ffff6a13ec3 in G4Radioactivation::DecayIt(G4Track const&, G4Step const&) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4processes.so
No symbol table info available.
#9  0x00007ffff3f99339 in G4SteppingManager::InvokePSDIP(unsigned long) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4tracking.so
No symbol table info available.
#10 0x00007ffff3f9960b in G4SteppingManager::InvokePostStepDoItProcs() () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4tracking.so
No symbol table info available.
#11 0x00007ffff3f99c76 in G4SteppingManager::Stepping() () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4tracking.so
No symbol table info available.
#12 0x00007ffff3fad6eb in G4TrackingManager::ProcessOneTrack(G4Track*) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4tracking.so
No symbol table info available.
#13 0x00007ffff75ab42d in G4EventManager::DoProcessing(G4Event*) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4event.so
No symbol table info available.
#14 0x00007ffff76cf9fc in G4WorkerRunManager::DoEventLoop(int, char const*, int) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4run.so
No symbol table info available.
#15 0x00007ffff767e112 in G4RunManager::BeamOn(int, char const*, int) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4run.so
No symbol table info available.
#16 0x00007ffff76d36b5 in G4WorkerRunManager::DoWork() () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4run.so
No symbol table info available.
#17 0x00007ffff7674a90 in G4MTRunManagerKernel::StartThread(G4WorkerThread*) () from /home/elias/software/geant4/geant4-v11.2.2-install/lib/libG4run.so
No symbol table info available.
#18 0x00007ffff4501b65 in std::execute_native_thread_routine (__p=<optimized out>) at ../../../../../libstdc++-v3/src/c++11/thread.cc:104
        __t = {_M_t = {<std::__uniq_ptr_impl<std::thread::_State, std::default_delete<std::thread::_State> >> = {
              _M_t = {<std::_Tuple_impl<0, std::thread::_State*, std::default_delete<std::thread::_State> >> = {<std::_Tuple_impl<1, std::default_delete<std::thread::_State> >> = {<std::_Head_base<1, std::def--Type <RET> for more, q to quit, c to continue without paging--
ault_delete<std::thread::_State>, true>> = {_M_head_impl = {<No data fields>}}, <No data fields>}, <std::_Head_base<0, std::thread::_State*, false>> = {
                    _M_head_impl = 0x555558016bd0}, <No data fields>}, <No data fields>}}, <No data fields>}}
#19 0x00007ffff4294a94 in start_thread (arg=<optimized out>) at ./nptl/pthread_create.c:447
        ret = <optimized out>
        pd = <optimized out>
        out = <optimized out>
        unwind_buf = {cancel_jmp_buf = {{jmp_buf = {140737194346176, -1922442972844585719, 140737194346176, -1123576, 2, 140737488329200, -1922442972823614199, -1922429223514075895}, mask_was_saved = 0}}, 
          priv = {pad = {0x0, 0x0, 0x0, 0x0}, data = {prev = 0x0, cleanup = 0x0, canceltype = 0}}}
        not_first_call = <optimized out>
#20 0x00007ffff4321c3c in clone3 () at ../sysdeps/unix/sysv/linux/x86_64/clone3.S:78
No locals.

I am struggling to take this information and find where the segfault occurs. Could someone give me a tip on how to proceed?

Thanks again.

5

You’re not going to like it… The segfault happened at frame #0 of the traceback:

But you don’t have the necessary debugging symbols in the build libraries to get any farther. You did your Geant4 installation in the usual recommended way, with the default CMAKE_BUILD_TYPE=Release. To get the detailed debugging information for the traceback, you’ll need to rebuild your Geant4 with CMAKE_BUILD_TYPE=RelWithDebInfo or Debug. You will probably want to rebuild your executable as well, with the same CMAKE_BUILD_TYPE setting.

I’d also suggest running with a single thread (/run/numberOfThreads 1) in your macro file. It makes debugging much easier.

6

I have rebuilt with CMAKE_BUILD_TYPE=RelWithDebInfo. Using the same gdb method as in my last post (on 1 thread now) results in:

G4PolarizationTransition::SampleGammaTransition (this=this@entry=0x7fffec1b22d8, nucpol=nucpol@entry=0x7fffec304590, twoJ1=twoJ1@entry=-4, twoJ2=twoJ2@entry=-6, L0=<optimized out>, Lp=<optimized out>, 
    mpRatio=mpRatio@entry=0, cosTheta=@0x7ffff09089c8: 0.9101682027022342, phi=@0x7ffff09089d0: 3.0797649869619117)
    at /home/elias/software/geant4/geant4-v11.2.2/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4PolarizationTransition.cc:361
361         for(std::size_t kappa2=0; kappa2<(newPol[k2]).size(); ++kappa2) {

Then: (gdb) backtrace full

#0  G4PolarizationTransition::SampleGammaTransition (this=this@entry=0x7fffec1b22d8, nucpol=nucpol@entry=0x7fffec304590, twoJ1=twoJ1@entry=-4, twoJ2=twoJ2@entry=-6, L0=<optimized out>, 
    Lp=<optimized out>, mpRatio=mpRatio@entry=0, cosTheta=@0x7ffff09089c8: 0.9101682027022342, phi=@0x7ffff09089d0: 3.0797649869619117)
    at /home/elias/software/geant4/geant4-v11.2.2/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4PolarizationTransition.cc:361
        kappa2 = 0
        lastNonZero = -1
        k2 = 0
        pol = @0x7fffec3045a0: {<std::_Vector_base<std::vector<std::complex<double>, std::allocator<std::complex<double> > >, std::allocator<std::vector<std::complex<double>, std::allocator<std::complex<double> > > > >> = {
            _M_impl = {<std::allocator<std::vector<std::complex<double>, std::allocator<std::complex<double> > > >> = {<std::__new_allocator<std::vector<std::complex<double>, std::allocator<std::complex<double> > > >> = {<No data fields>}, <No data fields>}, <std::_Vector_base<std::vector<std::complex<double>, std::allocator<std::complex<double> > >, std::allocator<std::vector<std::complex<double>, std::allocator<std::complex<double> > > > >::_Vector_impl_data> = {_M_start = 0x7fffec2c6bf0, _M_finish = 0x7fffec2c6c08, _M_end_of_storage = 0x7fffec2c6c08}, <No data fields>}}, <No data fields>}
        newlength = 18446744073709551611
        newPol = {<std::_Vector_base<std::vector<std::complex<double>, std::allocator<std::complex<double> > >, std::allocator<std::vector<std::complex<double>, std::allocator<std::complex<double> > > > >> = {
            _M_impl = {<std::allocator<std::vector<std::complex<double>, std::allocator<std::complex<double> > > >> = {<std::__new_allocator<std::vector<std::complex<double>, std::allocator<std::complex<double> > > >> = {<No data fields>}, <No data fields>}, <std::_Vector_base<std::vector<std::complex<double>, std::allocator<std::complex<double> > >, std::allocator<std::vector<std::complex<double>, std::allocator<std::complex<double> > > > >::_Vector_impl_data> = {_M_start = 0x0, _M_finish = 0x0, _M_end_of_storage = 0x0}, <No data fields>}}, <No data fields>}
        cachePtr = 0x0
        lastNonEmptyK2 = <optimized out>
#1  0x00007ffff66a99b4 in G4GammaTransition::SampleDirection (this=0x7fffec1b22b0, nuc=<optimized out>, ratio=<optimized out>, twoJ1=-4, twoJ2=-6, mp=4)
    at /home/elias/software/geant4/geant4-v11.2.2/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4GammaTransition.cc:185
        mpRatio = 0
        L0 = <optimized out>
        Lp = <optimized out>
        cosTheta = 0.9101682027022342
        phi = 3.0797649869619117
        np = 0x7fffec304590
        sinTheta = -nan(0xffffafffffffc)
#2  0x00007ffff66a976f in G4GammaTransition::SampleTransition (this=0x7fffec1b22b0, nucleus=0x7ffff0908cb0, newExcEnergy=1.2747660000000001, mpRatio=0, JP1=-4, JP2=-6, MP=4, shell=<optimized out>, 
    isDiscrete=true, isGamma=true) at /home/elias/software/geant4/geant4-v11.2.2/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4GammaTransition.cc:105
        result = 0x0
        bond_energy = 0
        etrans = <optimized out>
        lv = {pp = {static ToleranceTicks = 100, data = {0, 0, 0}}, ee = 199294.92728338746}
        mass = 199294.20740938745
        part = 0x555555a5e680
        emass = <optimized out>
        ecm = <optimized out>
        bst = {static ToleranceTicks = 100, data = {6.9533430133877208e-310, 0.2000018611602172, 8.654284465188766e-156}}
        energy = <optimized out>
        mom = <optimized out>
        res4mom = {pp = {static ToleranceTicks = 100, data = {6.9533430133782348e-310, 6.9533478658336494e-310, 0}}, ee = 0}
#3  0x00007ffff66ac78b in G4PhotonEvaporation::GenerateGamma (this=this@entry=0x7fffec1b21c0, nucleus=nucleus@entry=0x7ffff0908cb0)
    at /home/elias/software/geant4/geant4-v11.2.2/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4PhotonEvaporation.cc:543
        result = 0x0
        eexc = <optimized out>
        time = <optimized out>
        efinal = <optimized out>
        ratio = <optimized out>
        JP1 = <optimized out>
        JP2 = <optimized out>
        multiP = <optimized out>
        isGamma = <optimized out>
        isDiscrete = <optimized out>
--Type <RET> for more, q to quit, c to continue without paging--
        level = <optimized out>
        ntrans = <optimized out>
        isLL = <optimized out>
#4  0x00007ffff66adf2e in G4PhotonEvaporation::EmittedFragment (this=0x7fffec1b21c0, nucleus=0x7ffff0908cb0)
    at /home/elias/software/geant4/geant4-v11.2.2/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4PhotonEvaporation.cc:145
        fNucPStore = 0x7fffec2c5a50
        gamma = <optimized out>
#5  0x00007ffff6a75401 in G4ITDecay::DecayIt (this=0x7fffec1b2370) at /home/elias/software/geant4/geant4-v11.2.2/source/processes/hadronic/models/radioactive_decay/src/G4ITDecay.cc:87
        atRest = {pp = {static ToleranceTicks = 100, data = {0, 0, 0}}, ee = 199294.92728338746}
        parentParticle = {theMomentumDirection = {static ToleranceTicks = 100, data = {1, 0, 0}}, thePolarization = {static ToleranceTicks = 100, data = {0, 0, 0}}, 
          theParticleDefinition = 0x7fffec2d0f50, theElectronOccupancy = 0x0, thePreAssignedDecayProducts = 0x0, primaryParticle = 0x0, theKineticEnergy = 0, 
          theLogKineticEnergy = 1.7976931348623157e+308, theBeta = -1, theProperTime = 0, theDynamicalMass = 199294.92728338746, theDynamicalCharge = 84, theDynamicalSpin = 0, 
          theDynamicalMagneticMoment = 0, thePreAssignedDecayTime = -1, verboseLevel = 1, thePDGcode = 0}
        products = 0x7fffec304a90
        parentNucleus = {theA = 214, theZ = 84, theL = 0, theExcitationEnergy = 1.9946399999898858, theGroundStateMass = 199292.93264338747, theMomentum = {pp = {static ToleranceTicks = 100, data = {0, 
                0, 0}}, ee = 199294.92728338746}, thePolarization = 0x7fffec304590, creatorModel = -1, numberOfParticles = 0, numberOfCharged = 0, numberOfHoles = 0, numberOfChargedHoles = 0, 
          numberOfShellElectrons = 0, xLevel = 0, theParticleDefinition = 0x0, spin = 0, theCreationTime = 0, isLongLived = false}
        eOrGamma = <optimized out>
        theIonTable = <optimized out>
        daughterIon = <optimized out>
        dynDaughter = <optimized out>
#6  0x00007ffff6a7b58a in G4RadioactiveDecay::DoDecay (this=this@entry=0x7fffec1ac550, theParticleDef=..., theDecayTable=<optimized out>)
    at /home/elias/software/geant4/geant4-v11.2.2/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc:1138
        products = 0x0
        parentPlusQ = <optimized out>
        theDecayChannel = 0x7fffec268210
#7  0x00007ffff6a7bbce in G4RadioactiveDecay::DecayAnalog (this=0x7fffec1ac550, theTrack=..., decayTable=<optimized out>)
    at /home/elias/software/geant4/geant4-v11.2.2/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc:990
        theParticle = 0x7fffec2fd7c0
        theParticleDef = 0x7fffec2d0f50
        products = <optimized out>
        energyDeposit = <optimized out>
        finalGlobalTime = <optimized out>
        finalLocalTime = <optimized out>
        ParentEnergy = <optimized out>
        ParentDirection = {static ToleranceTicks = 100, data = {4.6355709537047487e-310, 6.9533393164241358e-310, 6.9533430134493802e-310}}
        numberOfSecondaries = <optimized out>
        modelID_forIT = <optimized out>
        modelID = <optimized out>
        modelID_forAtomicRelaxation = <optimized out>
#8  0x00007ffff6a89a36 in G4Radioactivation::DecayIt (this=0x7fffec1ac550, theTrack=...)
    at /home/elias/software/geant4/geant4-v11.2.2/source/processes/hadronic/models/radioactive_decay/src/G4Radioactivation.cc:871
        theParticle = <optimized out>
        theParticleDef = 0x7fffec2d0f50
        theDecayTable = 0x7fffec2fef30
#9  0x00007ffff420df72 in G4SteppingManager::InvokePSDIP (this=this@entry=0x7fffec04c7b0, np=np@entry=2) at /home/elias/software/geant4/geant4-v11.2.2/source/tracking/src/G4SteppingManager.cc:818
No locals.
#10 0x00007ffff420e08b in G4SteppingManager::InvokePostStepDoItProcs (this=this@entry=0x7fffec04c7b0) at /home/elias/software/geant4/geant4-v11.2.2/source/tracking/src/G4SteppingManager.cc:790
        Cond = <optimized out>
        np = 2
#11 0x00007ffff420e536 in G4SteppingManager::Stepping (this=0x7fffec04c7b0) at /home/elias/software/geant4/geant4-v11.2.2/source/tracking/src/G4SteppingManager.cc:225
        GeomStepLength = 0
        regionalAction = <optimized out>
#12 0x00007ffff42222c0 in G4TrackingManager::ProcessOneTrack (this=0x7fffec04b240, apValueG4Track=apValueG4Track@entry=0x7fffec300710)
--Type <RET> for more, q to quit, c to continue without paging--
    at /home/elias/software/geant4/geant4-v11.2.2/source/tracking/src/G4TrackingManager.cc:133
No locals.
#13 0x00007ffff761c5a5 in G4EventManager::DoProcessing (this=0x7fffec018080, anEvent=<optimized out>) at /home/elias/software/geant4/geant4-v11.2.2/source/event/src/G4EventManager.cc:211
        aTrajectory = <optimized out>
        secondaries = <optimized out>
        partDef = <optimized out>
        particleTrackingManager = <optimized out>
        previousTrajectory = 0x0
        currentState = <optimized out>
        center = {static ToleranceTicks = 100, data = {0, 0, 0}}
        navigator = <optimized out>
        track = 0x7fffec300710
        istop = <optimized out>
        trackingManagersToFlush = {
          _M_h = {<std::__detail::_Hashtable_base<G4VTrackingManager*, G4VTrackingManager*, std::__detail::_Identity, std::equal_to<G4VTrackingManager*>, std::hash<G4VTrackingManager*>, std::__detail::_Mod_range_hashing, std::__detail::_Default_ranged_hash, std::__detail::_Hashtable_traits<false, true, true> >> = {<std::__detail::_Hash_code_base<G4VTrackingManager*, G4VTrackingManager*, std::__detail::_Identity, std::hash<G4VTrackingManager*>, std::__detail::_Mod_range_hashing, std::__detail::_Default_ranged_hash, false>> = {<std::__detail::_Hashtable_ebo_helper<1, std::hash<G4VTrackingManager*>, true>> = {<std::hash<G4VTrackingManager*>> = {<std::__hash_base<unsigned long, G4VTrackingManager*>> = {<No data fields>}, <No data fields>}, <No data fields>}, <No data fields>}, <std::__detail::_Hashtable_ebo_helper<0, std::equal_to<G4VTrackingManager*>, true>> = {<std::equal_to<G4VTrackingManager*>> = {<std::binary_function<G4VTrackingManager*, G4VTrackingManager*, bool>> = {<No data fields>}, <No data fields>}, <No data fields>}, <No data fields>}, <std::__detail::_Map_base<G4VTrackingManager*, G4VTrackingManager*, std::allocator<G4VTrackingManager*>, std::__detail::_Identity, std::equal_to<G4VTrackingManager*>, std::hash<G4VTrackingManager*>, std::__detail::_Mod_range_hashing, std::__detail::_Default_ranged_hash, std::__detail::_Prime_rehash_policy, std::__detail::_Hashtable_traits<false, true, true>, true>> = {<No data fields>}, <std::__detail::_Insert<G4VTrackingManager*, G4VTrackingManager*, std::allocator<G4VTrackingManager*>, std::__detail::_Identity, std::equal_to<G4VTrackingManager*>, std::hash<G4VTrackingManager*>, std::__detail::_Mod_range_hashing, std::__detail::_Default_ranged_hash, std::__detail::_Prime_rehash_policy, std::__detail::_Hashtable_traits<false, true, true>, true>> = {<std::__detail::_Insert_base<G4VTrackingManager*, G4VTrackingManager*, std::allocator<G4VTrackingManager*>, std::__detail::_Identity, std::equal_to<G4VTrackingManager*>, std::hash<G4VTrackingManager*>, std::__detail::_Mod_range_hashing, std::__detail::_Default_ranged_hash, std::__detail::_Prime_rehash_policy, std::__detail::_Hashtable_traits<false, true, true> >> = {<No data fields>}, <No data fields>}, <std::__detail::_Rehash_base<G4VTrackingManager*, G4VTrackingManager*, std::allocator<G4VTrackingManager*>, std::__detail::_Identity, std::equal_to<G4VTrackingManager*>, std::hash<G4VTrackingManager*>, std::__detail::_Mod_range_hashing, std::__detail::_Default_ranged_hash, std::__detail::_Prime_rehash_policy, std::__detail::_Hashtable_traits<false, true, true>, std::integral_constant<bool, true> >> = {<No data fields>}, <std::__detail::_Equality<G4VTrackingManager*, G4VTrackingManager*, std::allocator<G4VTrackingManager*>, std::__detail::_Identity, std::equal_to<G4VTrackingManager*>, std::hash<G4VTrackingManager*>, std::__detail::_Mod_range_hashing, std::__detail::_Default_ranged_hash, std::__detail::_Prime_rehash_policy, std::__detail::_Hashtable_traits<false, true, true>, true>> = {<No data fields>}, <std::__detail::_Hashtable_alloc<std::allocator<std::__detail::_Hash_node<G4VTrackingManager*, false> > >> = {<std::__detail::_Hashtable_ebo_helper<0, std::allocator<std::__detail::_Hash_node<G4VTrackingManager*, false> >, true>> = {<std::allocator<std::__detail::_Hash_node<G4VTrackingManager*, false> >> = {<std::__new_allocator<std::__detail::_Hash_node<G4VTrackingManager*, false> >> = {<No data fields>}, <No data fields>}, <No data fields>}, <No data fields>}, <std::_Enable_default_constructor<true, std::__detail::_Hash_node_base>> = {<No data fields>}, _M_buckets = 0x7ffff09093f0, _M_bucket_count = 1, 
            _M_before_begin = {_M_nxt = 0x0}, _M_element_count = 0, _M_rehash_policy = {static _S_growth_factor = 2, _M_max_load_factor = 1, _M_next_resize = 0}, _M_single_bucket = 0x0}}
        nses = <optimized out>
#14 0x00007ffff761d190 in G4EventManager::ProcessOneEvent (this=<optimized out>, anEvent=<optimized out>) at /home/elias/software/geant4/geant4-v11.2.2/source/event/src/G4EventManager.cc:447
No locals.
#15 0x00007ffff772d3c4 in G4WorkerRunManager::ProcessOneEvent (this=0x7fffec04b000, i_event=<optimized out>) at /home/elias/software/geant4/geant4-v11.2.2/source/run/src/G4WorkerRunManager.cc:287
No locals.
#16 0x00007ffff772d5fe in G4WorkerRunManager::DoEventLoop (this=0x7fffec04b000, n_event=10000, macroFile=0x0, n_select=-1)
    at /home/elias/software/geant4/geant4-v11.2.2/source/run/src/G4WorkerRunManager.cc:270
        i_event = -1
#17 0x00007ffff76e0d72 in G4RunManager::BeamOn (this=0x7fffec04b000, n_event=10000, macroFile=0x0, n_select=-1) at /home/elias/software/geant4/geant4-v11.2.2/source/run/src/G4RunManager.cc:261
        cond = <optimized out>
#18 0x00007ffff7730dad in G4WorkerRunManager::DoWork (this=<optimized out>) at /home/elias/software/geant4/geant4-v11.2.2/source/run/src/G4WorkerRunManager.cc:721
        numevents = 10000
        skipInitialization = false
        cmds = {<std::_Vector_base<G4String, std::allocator<G4String> >> = {
            _M_impl = {<std::allocator<G4String>> = {<std::__new_allocator<G4String>> = {<No data fields>}, <No data fields>}, <std::_Vector_base<G4String, std::allocator<G4String> >::_Vector_impl_data> = {_M_start = 0x7fffec283290, _M_finish = 0x7fffec283330, _M_end_of_storage = 0x7fffec283330}, <No data fields>}}, <No data fields>}
        uimgr = <optimized out>
        macroFile = {<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >> = {static npos = 18446744073709551615, 
            _M_dataplus = {<std::allocator<char>> = {<std::__new_allocator<char>> = {<No data fields>}, <No data fields>}, _M_p = 0x7ffff09096c0 ""}, _M_string_length = 0, {
              _M_local_buf = "\000\357\b\354\377\177\000\000|R\272\367\377\177\000", _M_allocated_capacity = 140737153396480}}, <No data fields>}
        numSelect = <optimized out>
--Type <RET> for more, q to quit, c to continue without paging--
        mrm = 0x5555555d3590
        nextAction = <optimized out>
#19 0x00007ffff76d7c80 in G4MTRunManagerKernel::StartThread (context=<optimized out>) at /home/elias/software/geant4/geant4-v11.2.2/source/run/src/G4MTRunManagerKernel.cc:210
        masterRM = 0x5555555d3590
        thisID = <optimized out>
        masterEngine = <optimized out>
        wrm = 0x7fffec04b000
        wrmm = {<std::unique_lock<std::mutex>> = {_M_device = 0x7ffff77605c0 <(anonymous namespace)::workerRMMutex>, _M_owns = false}, <No data fields>}
        detector = <optimized out>
        physicslist = <optimized out>
#20 0x00007ffff4775b65 in std::execute_native_thread_routine (__p=<optimized out>) at ../../../../../libstdc++-v3/src/c++11/thread.cc:104
        __t = {_M_t = {<std::__uniq_ptr_impl<std::thread::_State, std::default_delete<std::thread::_State> >> = {
              _M_t = {<std::_Tuple_impl<0, std::thread::_State*, std::default_delete<std::thread::_State> >> = {<std::_Tuple_impl<1, std::default_delete<std::thread::_State> >> = {<std::_Head_base<1, std::default_delete<std::thread::_State>, true>> = {_M_head_impl = {<No data fields>}}, <No data fields>}, <std::_Head_base<0, std::thread::_State*, false>> = {
                    _M_head_impl = 0x555558016e10}, <No data fields>}, <No data fields>}}, <No data fields>}}
#21 0x00007ffff4508a94 in start_thread (arg=<optimized out>) at ./nptl/pthread_create.c:447
        ret = <optimized out>
        pd = <optimized out>
        out = <optimized out>
        unwind_buf = {cancel_jmp_buf = {{jmp_buf = {140737230522048, -5570738301346945212, 140737230522048, -1123640, 2, 140737488327872, -5570738301359528124, -5570730055350917308}, 
              mask_was_saved = 0}}, priv = {pad = {0x0, 0x0, 0x0, 0x0}, data = {prev = 0x0, cleanup = 0x0, canceltype = 0}}}
        not_first_call = <optimized out>
#22 0x00007ffff4595c3c in clone3 () at ../sysdeps/unix/sysv/linux/x86_64/clone3.S:78
No locals.

I am again struggling to find the important information in this so I would be very thankful for some more input. Thanks for your help.

7

The traceback is the start of your debugging adventure, not the final answer. The segfault happened in “frame #0.” In your case, that’s at line 361 of G4PolarizationTransition.cc. In the debugger, you can see the source code around that failure point with list.

That said, there might be a clue in frame #1, it’s showing that one of the local variable, sinTheta = -nan. That doesn’t sound good.

From your original (non-debug) traceback, frame #6 has the specific particle being decayed. And I see in frame #5 theres a variable parentNucleus which tells you that it’s A=214, Z=84 (Po-214). I wonder if you’d be able to reproduce this by directly specifying Po-214 as your primary. That would certainly make debugging easier.

8

I think I may have found the problem now and it seems to be a bug in Geant4’s code. I will do my best to summarize from the beginning what I think happens.

  1. I am simulating nuclear decays where gammas are emitted. It is important for me that the gammas are correlated and that the angular distributions of gamma-ray cascades are accurate. Therefore, in my PhysicsList.cc, I am using:
    G4DeexPrecoParameters* deex = G4NuclearLevelData::GetInstance()->GetParameters();
    deex->SetCorrelatedGamma(true);
  1. Geant4 simulates a transition from an excited nuclear state to a lower one or a ground state. The states each have an angular momentum (J) and parity (P). The data for this is stored in YourInstallFolder\share\Geant4\data\PhotonEvaporation5.7. Here J and P information is stored as J*P. For example, J=2, P=+ is stored as 2 and J=3/2, P=- is stored as -1.5.
  2. The problem arises when the second nuclear state has negative parity. So P=- and the value for J*P is <0. Let us call this JP2 < 0.
  3. Geant4 uses Geant4SourceFolder\source\processes\hadronic\models\de_excitation\photon_evaporation\src\ to handle correlated gammas. Somewhere here JP2 is multiplied by two.
  4. In G4PolarizationTransition.cc, our 2*JP2 is stored in variable fTwoJ2. In the case above, it is negative. With verbosity, it is printed and valued as expected.
  5. This happens on line 279 in G4PolarizationTransition.cc:
  std::size_t newlength = fTwoJ2+1;
  1. But std::size_t is unsigned and fTwoJ2+1 can be negative. In the problematic case newlength caps out at ~2^64. This is seen in frame #0 in my last post where newlength=18446744073709551611.
  2. On line 287 in G4PolarizationTransition.cc there is a for loop over newlength, which now is very big:
for(G4int k2=0; k2<(G4int)newlength; ++k2) {
  1. From here on, bad stuff happens which I think causes the segmentation fault.

I have tested several examples and this segfault happens every time fTwoJ2+1<0. Also fTwoJ2+1=0 causes a segfault. So this segfault always happens when the parity of the second state is negative. It never happens when the parity of the second state is positive or J=0.

There are many nuclear decays where this happens.

It feels strange that this has gone undetected and I am not excluding the possibility that I am wrong. I am trying to learn both Geant4 and C++ here, so please inform me of any mistakes.

9

I think you’ve made that really clear, Elias. @ribon is our Hadronics WG coordinator, @civanch is our EM expert, including the de-excitation stuff. This looks like a candidate for a Bugzilla report, but I’d like them to chime in first, in case I’m missing something.

1 Like
10

In the meantime, I did a fix that appears to remove these segmentation faults.

In geant4-v11.3.0\source\processes\hadronic\models\de_excitation\photon_evaporation\src\G4PolarizationTransition.cc I changed lines 236-237 from

  fTwoJ1 = twoJ1;  // add abs to remove negative J
  fTwoJ2 = twoJ2;

to

  fTwoJ1 = std::abs(twoJ1);
  fTwoJ2 = std::abs(twoJ2);

Then I rebuilt Geant4 and it works.

Curiously, there was already an existing comment about adding an absolute value.

And, of course, a big thank you Michael for your help.

11

Would you be willing to create a bug report, including your fix above? I am especially amused by the comment :-/

12

Yep, I will fill out a bug report

13

I have made a bug report now and marked my fix as the solution in this thread :+1:

14

Dear @Elias_Arnqvist @mkelsey

Thank you for reporting this bug!

Out of curiosity, does ignoring parity have any impact further down the code, such as on the generation of the spatial distribution or selection rules?

Perhaps Geant4 developers can propose an extended patch which adapt the code to handle negative parities rather than simply ignoring them.

Thank you for your time.

Best,
Alvaro

15

Hi @atolosad

I was also initially worried about removing parity information, but it seems the parity is not used past that point. It is as far down the code as the parity can be ignored while removing the segmentation faults. The angular momentum is used past this point to determine the angular distribution of gammas, but I do not think parity is used for this. The naming convention also switches from JP2 to twoJ2.

I agree that it would be best for the Geant4 developers to take a look at what the best way to do this is :slight_smile:

16

This topic was automatically closed 2 days after the last reply. New replies are no longer allowed.