Simple Gamma Detection Simulation

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: 10.7.4
Operating System: RHEL9
Compiler/Version: 11.4.1
CMake Version: 3.20.2

I’ve simulated a Cs137 point source sitting next to a NaI crystal and recorded the energy deposited in the crystal to form an energy spectra. I’ve uploaded the results of the simulated energy spectra along with one collected from a NaI detector I have in my lab. The big noticable difference is in the Compton back scatter peak around 200keV which is absent from the simulation. Can someone help pin point how to figure out if I’m using the right physics libraries to properly simulate an energy spectra of a gamma less than 2MeV?

The simulation was done using the topasmc.org geant4 framework.

Here is a cut and paste of the simulation output at the initialization stage listing all the physics processes.

Registered Physics Processes:
     Transportation,        StepLimiter,     UserSpecialCut,               phot
              compt,               conv,               Rayl,                msc
              eIoni,              eBrem,        CoulombScat,                msc
              eIoni,              eBrem,            annihil,        CoulombScat
                msc,            ionIoni,                msc,             muIoni
        CoulombScat,             muIoni,                msc,              hIoni
        CoulombScat,              hIoni,                msc,              hIoni
        CoulombScat,              hIoni,                msc,              hIoni
        CoulombScat,                msc,              hIoni,        CoulombScat
              hIoni,              hIoni,                msc,            ionIoni
                msc,            ionIoni,RadioactiveDecayBase

Visualization verbosity changed to quiet (0)
My Begin Run has been called edited

### ===  Deexcitation model UAtomDeexcitation is activated for 1 region:
          DefaultRegionForTheWorld  1  1  0
### ===  Auger cascade flag: 1
### ===  Ignore cuts flag:   1

phot:  for gamma SubType=12 BuildTable=0
      LambdaPrime table from 200 keV to 500 MeV in 25 bins 
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
 LivermorePhElectric : Emin=    0 eV  Emax=  500 MeV  SauterGavrila Fluo

compt:  for gamma SubType=13 BuildTable=1
      Lambda table from 100 eV  to 1 MeV, 7 bins/decade, spline: 1
      LambdaPrime table from 1 MeV to 500 MeV in 20 bins 
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
       Klein-Nishina : Emin=    0 eV  Emax=  500 MeV

conv:  for gamma SubType=14 BuildTable=1
      Lambda table from 1.022 MeV to 500 MeV, 33 bins/decade, spline: 1
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
     BetheHeitlerLPM : Emin=    0 eV  Emax=  500 MeV  ModifiedTsai

Rayl:  for gamma SubType=11 BuildTable=1
      Lambda table from 100 eV  to 100 keV, 7 bins/decade, spline: 0
      LambdaPrime table from 100 keV to 500 MeV in 27 bins 
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
   LivermoreRayleigh : Emin=    0 eV  Emax=  500 MeV  CullenGenerator

msc:  for e-  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            UrbanMsc : Emin=    0 eV  Emax=  100 MeV Nbins=42 100 eV  - 100 MeV
              StepLim=UseSafety Rfact=0.04 Gfact=2.5 Sfact=0.6 DispFlag:1 Skin=1 Llimit=1
        WentzelVIUni : Emin=  100 MeV Emax=  500 MeV Nbins=7 100 MeV - 500 MeV
              StepLim=UseSafety Rfact=0.04 Gfact=2.5 Sfact=0.6 DispFlag:1 Skin=1 Llimit=1

eIoni:  for e-  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 1 mm), integ: 1, fluct: 1, linLossLim= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        MollerBhabha : Emin=    0 eV  Emax=  500 MeV

eBrem:  for e-  SubType=3
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      LPM flag: 1 for E > 0.5 GeV,  VertexHighEnergyTh(GeV)= 100000
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
             eBremSB : Emin=    0 eV  Emax=  500 MeV  ModifiedTsai

CoulombScat:  for e-, integral:1  SubType=1 BuildTable=1
      Lambda table from 100 MeV to 500 MeV, 7 bins/decade, spline: 1
      ThetaMin(p) < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering : Emin=  100 MeV Emax=  500 MeV

msc:  for e+  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            UrbanMsc : Emin=    0 eV  Emax=  100 MeV Nbins=42 100 eV  - 100 MeV
              StepLim=UseSafety Rfact=0.04 Gfact=2.5 Sfact=0.6 DispFlag:1 Skin=1 Llimit=1
        WentzelVIUni : Emin=  100 MeV Emax=  500 MeV Nbins=7 100 MeV - 500 MeV
              StepLim=UseSafety Rfact=0.04 Gfact=2.5 Sfact=0.6 DispFlag:1 Skin=1 Llimit=1

eIoni:  for e+  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 1 mm), integ: 1, fluct: 1, linLossLim= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        MollerBhabha : Emin=    0 eV  Emax=  500 MeV

eBrem:  for e+  SubType=3
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      LPM flag: 1 for E > 0.5 GeV,  VertexHighEnergyTh(GeV)= 100000
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
             eBremSB : Emin=    0 eV  Emax=  500 MeV  ModifiedTsai

annihil:  for e+, integral:1  SubType=5 BuildTable=0
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            eplus2gg : Emin=    0 eV  Emax=  500 MeV

CoulombScat:  for e+, integral:1  SubType=1 BuildTable=1
      Lambda table from 100 MeV to 500 MeV, 7 bins/decade, spline: 1
      ThetaMin(p) < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering : Emin=  100 MeV Emax=  500 MeV

msc:  for proton  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        WentzelVIUni : Emin=    0 eV  Emax=  500 MeV Nbins=49 100 eV  - 500 MeV
              StepLim=Minimal Rfact=0.2 Gfact=2.5 Sfact=0.6 DispFlag:0 Skin=1 Llimit=1

hIoni:  for proton  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 0.1 mm), integ: 1, fluct: 1, linLossLim= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
               Bragg : Emin=    0 eV  Emax=    2 MeV
          BetheBloch : Emin=    2 MeV Emax=  500 MeV

CoulombScat:  for proton, integral:1  SubType=1 BuildTable=1
      Lambda table from threshold  to 500 MeV, 7 bins/decade, spline: 1
      ThetaMin(p) < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering : Emin=    0 eV  Emax=  500 MeV

msc:  for GenericIon  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            UrbanMsc : Emin=    0 eV  Emax=  500 MeV
              StepLim=Minimal Rfact=0.2 Gfact=2.5 Sfact=0.6 DispFlag:0 Skin=1 Llimit=1

ionIoni:  for GenericIon  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 0.1 mm), integ: 1, fluct: 1, linLossLim= 0.02
      Stopping Power data for 17 ion/material pairs
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            BraggIon : Emin=    0 eV  Emax=    2 MeV
          BetheBloch : Emin=    2 MeV Emax=  500 MeV
======================================================================
======          Radioactive Decay Physics Parameters           =======
======================================================================
Max life time                                     1.4427e+06 ps
Internal e- conversion flag                       1
Stored internal conversion coefficients           1
Enable correlated gamma emission                  0
Max 2J for sampling of angular correlations       10
Atomic de-excitation enabled                      1
Auger electron emission enabled                   1
Auger cascade enabled                             1
Check EM cuts disabled for atomic de-excitation   1
Use Bearden atomic level energies                 0
======================================================================

msc:  for alpha  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            UrbanMsc : Emin=    0 eV  Emax=  500 MeV Nbins=49 100 eV  - 500 MeV
              StepLim=Minimal Rfact=0.2 Gfact=2.5 Sfact=0.6 DispFlag:0 Skin=1 Llimit=1

ionIoni:  for alpha  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 0.1 mm), integ: 1, fluct: 1, linLossLim= 0.02
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            BraggIon : Emin=    0 eV  Emax=7.9452 MeV
          BetheBloch : Emin=7.9452 MeV Emax=  500 MeV

msc:  for anti_proton  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        WentzelVIUni : Emin=    0 eV  Emax=  500 MeV Nbins=49 100 eV  - 500 MeV
              StepLim=Minimal Rfact=0.2 Gfact=2.5 Sfact=0.6 DispFlag:0 Skin=1 Llimit=1

hIoni:  for anti_proton  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 0.1 mm), integ: 1, fluct: 1, linLossLim= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            ICRU73QO : Emin=    0 eV  Emax=    2 MeV
          BetheBloch : Emin=    2 MeV Emax=  500 MeV

CoulombScat:  for anti_proton, integral:1  SubType=1 BuildTable=1
      Lambda table from threshold  to 500 MeV, 7 bins/decade, spline: 1
      ThetaMin(p) < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering : Emin=    0 eV  Emax=  500 MeV

msc:  for kaon+  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        WentzelVIUni : Emin=    0 eV  Emax=  500 MeV Nbins=49 100 eV  - 500 MeV
              StepLim=Minimal Rfact=0.2 Gfact=2.5 Sfact=0.6 DispFlag:0 Skin=1 Llimit=1

hIoni:  for kaon+  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 0.1 mm), integ: 1, fluct: 1, linLossLim= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
               Bragg : Emin=    0 eV  Emax=1.05231 MeV
          BetheBloch : Emin=1.05231 MeV Emax=  500 MeV

CoulombScat:  for kaon+, integral:1  SubType=1 BuildTable=1
      Lambda table from threshold  to 500 MeV, 7 bins/decade, spline: 1
      ThetaMin(p) < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering : Emin=    0 eV  Emax=  500 MeV

msc:  for kaon-  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        WentzelVIUni : Emin=    0 eV  Emax=  500 MeV Nbins=49 100 eV  - 500 MeV
              StepLim=Minimal Rfact=0.2 Gfact=2.5 Sfact=0.6 DispFlag:0 Skin=1 Llimit=1

hIoni:  for kaon-  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 0.1 mm), integ: 1, fluct: 1, linLossLim= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            ICRU73QO : Emin=    0 eV  Emax=1.05231 MeV
          BetheBloch : Emin=1.05231 MeV Emax=  500 MeV

CoulombScat:  for kaon-, integral:1  SubType=1 BuildTable=1
      Used Lambda table of kaon+
      ThetaMin(p) < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering : Emin=    0 eV  Emax=  500 MeV

msc:  for mu+  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        WentzelVIUni : Emin=    0 eV  Emax=  500 MeV Nbins=49 100 eV  - 500 MeV
              StepLim=Minimal Rfact=0.2 Gfact=2.5 Sfact=0.6 DispFlag:0 Skin=1 Llimit=1

muIoni:  for mu+  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 0.1 mm), integ: 1, fluct: 1, linLossLim= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
               Bragg : Emin=    0 eV  Emax=  200 keV
          BetheBloch : Emin=  200 keV Emax=  500 MeV

CoulombScat:  for mu+, integral:1  SubType=1 BuildTable=1
      Lambda table from threshold  to 500 MeV, 7 bins/decade, spline: 1
      ThetaMin(p) < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering : Emin=    0 eV  Emax=  500 MeV

msc:  for mu-  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        WentzelVIUni : Emin=    0 eV  Emax=  500 MeV Nbins=49 100 eV  - 500 MeV
              StepLim=Minimal Rfact=0.2 Gfact=2.5 Sfact=0.6 DispFlag:0 Skin=1 Llimit=1

muIoni:  for mu-  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 0.1 mm), integ: 1, fluct: 1, linLossLim= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            ICRU73QO : Emin=    0 eV  Emax=  200 keV
          BetheBloch : Emin=  200 keV Emax=  500 MeV

CoulombScat:  for mu-, integral:1  SubType=1 BuildTable=1
      Used Lambda table of mu+
      ThetaMin(p) < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering : Emin=    0 eV  Emax=  500 MeV

msc:  for pi+  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        WentzelVIUni : Emin=    0 eV  Emax=  500 MeV Nbins=49 100 eV  - 500 MeV
              StepLim=Minimal Rfact=0.2 Gfact=2.5 Sfact=0.6 DispFlag:0 Skin=1 Llimit=1

hIoni:  for pi+  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 0.1 mm), integ: 1, fluct: 1, linLossLim= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
               Bragg : Emin=    0 eV  Emax=297.505 keV
          BetheBloch : Emin=297.505 keV Emax=  500 MeV

CoulombScat:  for pi+, integral:1  SubType=1 BuildTable=1
      Lambda table from threshold  to 500 MeV, 7 bins/decade, spline: 1
      ThetaMin(p) < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering : Emin=    0 eV  Emax=  500 MeV

msc:  for pi-  SubType= 10
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        WentzelVIUni : Emin=    0 eV  Emax=  500 MeV Nbins=49 100 eV  - 500 MeV
              StepLim=Minimal Rfact=0.2 Gfact=2.5 Sfact=0.6 DispFlag:0 Skin=1 Llimit=1

hIoni:  for pi-  SubType=2
      dE/dx and range tables from 100 eV  to 500 MeV in 49 bins
      Lambda tables from threshold to 500 MeV, 7 bins/decade, spline: 1
      StepFunction=(0.2, 0.1 mm), integ: 1, fluct: 1, linLossLim= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            ICRU73QO : Emin=    0 eV  Emax=297.505 keV
          BetheBloch : Emin=297.505 keV Emax=  500 MeV

CoulombScat:  for pi-, integral:1  SubType=1 BuildTable=1
      Used Lambda table of pi+
      ThetaMin(p) < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering : Emin=    0 eV  Emax=  500 MeV

=======================================================================
--> G4MTRunManager::CreateAndStartWorkers() --> Initializing workers...
=======================================================================

G4WT0 >      Transportation,        StepLimiter,     UserSpecialCut,               phot
G4WT0 >               compt,               conv,               Rayl,                msc
G4WT0 >               eIoni,              eBrem,        CoulombScat,                msc
G4WT0 >               eIoni,              eBrem,            annihil,        CoulombScat
G4WT0 >                 msc,            ionIoni,                msc,             muIoni
G4WT0 >         CoulombScat,             muIoni,                msc,              hIoni
G4WT0 >         CoulombScat,              hIoni,                msc,              hIoni
G4WT0 >         CoulombScat,              hIoni,                msc,              hIoni
G4WT0 >         CoulombScat,                msc,              hIoni,        CoulombScat
G4WT0 >               hIoni,              hIoni,                msc,            ionIoni
G4WT0 >                 msc,            ionIoni,RadioactiveDecayBase
G4WT0 > Begin processing for Run: 0, History: 0


The backscatter peak (normally a broad bump) is a real-world experimental effect caused by gamma rays scattering at large angles (very roughly 180 degrees) from materials surrounding the detector into the detector. Your simulation only has the detector and a source so I would not expect a backscatter peak.

1 Like

This is very useful. I’ll add the aluminum casing used to insulate the NaI crystal and surrounding lead shielding. Let’s see if the Compton back-scatter appears.