Simulating 133-Ba Decay

Hi All,

I am currently attempting to simulate the decay of 133-Ba. I am familiar with the the basic ways to generate particles. Would any of you happen to know how to setup a source based on the decay of 133-Ba?

Currently, I have attempted to use G4ParticleGun and its member functions to generate x-rays, with probabilities based off of tables I have found for the decay: https://www.nndc.bnl.gov/nudat3/decaysearchdirect.jsp?nuc=133Ba&unc=NDS

Thank you all in advance for your help.

There are two ways you could do it. If you want to let Geant4 handle all of the complications, you can use

/gun/ion 133 56

and Geant4 will do the radioactive decays, generating gammas with all the right probabilities.

If you want to do it yourself, you are probably better off using GPS (G4GeneralParticleSource) as your generator, rather than G4ParticleGun. You can set up energy spectra via histograms.

why not to start directly with Ba133 ?
Here, a macro to be run with rdecay01 (in examples/extended/radioactivedecay), to see the potentialities of Geant4.

Ba133.mac.txt (195 Bytes)

Thank you so much for your quick response.

When passing the command: /gun/ion 133 56 via a macro, is the ion placed at the origin? I currently have my method starting the source particle at the origin and traveling toward the detector geometry.

Also, do any other headers/libraries need to be included to use this command?

Thanks again!

Thank you for your quick response. I noticed in this example that there is a lot of other functionality included, is there a way to use the methods of this example without the use of the other functionality/recording/verbosity?

If you have your own primary action, and are just calling into your G4ParticleGun instance, you may not want to use the macro commands. Instead, you can instantiate the Ba-133 source in your code via:

G4ParticleDefinition* ion = G4IonTable::GetIonTable()->GetIon(133, 56, 0.); fParticleGun->SetParticleDefinition(ion);

where “fParticleGun” would be whatever you use for that data member.

If you use something like the code above, you’ll need to include the headers appropriate for it. If you switch to using macro commands, the macro command will already have the code defined and compiled for you.

Thank you again so much for your quick responses.

In my primary generator action, I currently have:

#include "PrimaryGeneratorAction.hh"

#include "G4LogicalVolumeStore.hh"
#include "G4LogicalVolume.hh"
#include "G4Box.hh"
#include "G4RunManager.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
#include "G4TransportationManager.hh" 

#include "EventHandler.hh"
#include "IO.hh"
#include "Storage.hh"

...
G4int n_particle = 1;
fParticleGun  = new G4ParticleGun(n_particle);

G4ParticleDefinition* particle ; 
particle  = particleTable->FindParticle(particleName="gamma");
energy = 31*keV;
position  = G4ThreeVector(0*mm,0.*mm,0.*mm) ; 
direction = G4ThreeVector(1,0,0) ;

fParticleGun->SetParticlePosition(position);
fParticleGun->SetParticleDefinition(particle);
fParticleGun->SetParticleMomentumDirection(direction);
fParticleGun->SetParticleEnergy(energy);
fParticleGun->SetParticleMomentumDirection(direction);
fParticleGun->GeneratePrimaryVertex(anEvent);
...

Would I then replace the above definitions as necessary given what you suggest:

#include "PrimaryGeneratorAction.hh"

#include "G4LogicalVolumeStore.hh"
#include "G4LogicalVolume.hh"
#include "G4Box.hh"
#include "G4RunManager.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
#include "G4TransportationManager.hh" 

#include "EventHandler.hh"
#include "IO.hh"
#include "Storage.hh"

...
G4int n_particle = 1;
fParticleGun  = new G4ParticleGun(n_particle);

G4ParticleDefinition* ion =  G4IonTable::GetIonTable()->GetIon(133, 56, 0.); 
position  = G4ThreeVector(0*mm,0.*mm,0.*mm) ; 


fParticleGun->SetParticlePosition(position);
fParticleGun->SetParticleDefinition(ion);
fParticleGun->GeneratePrimaryVertex(anEvent);
...

From the headers listed, do you know if these are sufficient, or are there others I should include? Please let me know if you need to know more information regarding their contents as some are home-brewed.

Thank you again so much for all of your help.

This is not the best practice, as it involves string comparisons which are quite inefficient. You can instantiate particles directly, because they are pure singletons; in your case, do

#include "G4Gamma.hh"
...
particle = G4Gamma::Definition();

In your replacement code, for the Ba-133 decays, the only additional thing you’ll need is

#include "G4IonTable.hh"

(which you can tell you need because you’re making use of G4IonTable functions in your code).

Thank you again for your help.

I tried implementing what we discussed and received the following error messages:

G4NucleiProperties::GetMassExccess: Wrong values for A = 56 and Z = 133
G4BinaryLightIonReaction no final state for: 
 Primary 0x7fe747b8ef90, (A,Z)=(56,133) , kinetic energy 1000
 Target nucleus (A,Z)=(1,1)
 if frequent, please submit above information as bug report


-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : GeomNav0003
      issued by : G4Navigator::ComputeStep()
Stuck Track: potential geometry or navigation problem.
  Track stuck, not moving for 25 steps.
  Current  phys volume: 'ReflectorLayer'
   - at position : (398.2,0,0)
     in direction: (1,0,0)
    (local position: (-3.061616997870903e-18,0,0.05000000000001137))
    (local direction: (6.123233995736766e-17,0,-1)).
  Previous phys volume: 'ReflectorLayer'

  Likely geometry overlap - else navigation problem !
 Track *abandoned* due to excessive number of Zero steps. Event aborted. 


*** Event Must Be Aborted ***
G4Track (0x7fe747b90240) - track ID = 1, parent ID = 0
 Particle type : E133-56 - creator process : not available
 Kinetic energy : 1 GeV - Momentum direction : (1,0,0)
 Step length : 0 fm  - total energy deposit : 0 eV 
 Pre-step point : (398.2,0,0) - Physical volume : ReflectorLayer (Polyvinyl Chloride)
 - defined by : ionInelastic - step status : 4
 Post-step point : (398.2,0,0) - Physical volume : ReflectorLayer (Polyvinyl Chloride)
 - defined by : ionInelastic - step status : 7
 *** Note: Step information might not be properly updated.

-------- EEEE -------- G4Exception-END --------- EEEE -------

Set the energy to zero. In your simulation, you want the Ba-133 to decay in place, and just emit the appropriate gamma rays. By setting the energy to 1 MeV, you’re basically throwing a beam of Ba-133 nuclei, and maybe don’t have the right physics list for that :slight_smile:

Thank you. I have set the energy to zero, but am still obtaining an error:

G4NucleiProperties::GetNuclearMass: Wrong values for A = 56 and Z = 133

Thank you for helping me to debug this. Do you think that I am not using the right physics list? Currently I am using G4EmLivermorePhysics and G4OpticalPhysics via:

G4VModularPhysicsList* physicsList = new FTFP_BERT;
physicsList->ReplacePhysics(new G4EmLivermorePhysics());
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
physicsList->RegisterPhysics(opticalPhysics);
runManager-> SetUserInitialization(physicsList);

Thanks again!

I think I gave you a bad bit of code. You should be able to debug it yourself :slight_smile: Look at the error message, which is saying you asked for a nucleus with A=56 (a total of 56 protons and neutrons combined) and Z=133 (133 protons).

Ah yes! I am sorry, I should have took the time to actually read the message. Was being lazy.

Thank you again for all of your help.

No worries! Hopefully it’s working for you now.

It compiles and executes fine, but I am not detecting any interaction with my scintillator material. Could it be that I need to set the energy spectra? Or is it just that the secondaries that are generated are not traveling through my geometry?

Where did you place the Ba-133 source? Radioactive decays are isotropic in Geant4. to get a decent event rate, you’ll need to have the source very close to your detector.

I placed the source 0.1 mm from the detector. I also checked in my visualization software and saw no production of any secondaries.

Are you sure that RadioactiveDecay process is registered in your PhysicsList ?

I don’t think it is. I have just tried in my main source file:

G4VModularPhysicsList* physicsList = new FTFP_BERT;
physicsList->ReplacePhysics(new G4EmLivermorePhysics());
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
G4RadioactiveDecay* decayPhysics = new G4RadioactiveDecay();
physicsList->RegisterPhysics(opticalPhysics);
physicsList->RegisterPhysics(decayPhysics);
runManager-> SetUserInitialization(physicsList);

But am running into an error:

error: cannot initialize a parameter of type 'G4VPhysicsConstructor *' with an lvalue of type 'G4RadioactiveDecay *'
physicsList->RegisterPhysics(decayPhysics);

I assume this is not the appropriate way to register the radioactive decay processes.

Nope :slight_smile: G4RadioactiveDecay is a process, not a “physics constructor.” What you want is

physicsList->RegisterPhysics(new G4RadioactiveDecayPhysics);

You don’t need to put the pointer into a local variable, unless you really want to. The physics list takes ownership, so you don’t have to hold onto the pointer for later deletion.

1 Like