Hi all. I’m attempting to create a simple germanium gamma detector script with an 800keV source as well as a Cs137 source. I’ve modelled the detector and created the sources closely following along with the example gamma-detector script, but am running into trouble getting accurate results out. I’d like to define the sources 30cm above the detector, but when I do the bins within ‘pulse-height’ tallies are blank or nearly blank. When I define the source inside of the crystal, it seems to work better, but the spectrum is still incorrect and shows peaks at the wrong energies.
I am using OpenMC version 0.14.0 installed using Docker and hosted on Jupyter Notebook.
Anyone have ideas of what I am doing wrong? I’ve looked into changing the domains of the sources with no success, I think it may have something to do with the tallies as defined, or maybe I fundamentally don’t understand how IndependentSource works. I would appreciate any input.
> ### IMPORTS ###
>
> import openmc
> import math
> from matplotlib import pyplot as plt
> import importlib
> import numpy as np
> import pandas as pd
> from ipywidgets import interact, widgets, fixed
> import openmc_source_plotter
> import warnings
> from scipy.constants import Avogadro
> warnings.filterwarnings('ignore')
> import os
> WorkingFolder = os.getcwd()
>
> ### DEF MATERIALS ###
>
> aluminum_pure = openmc.Material(0, 'pure aluminium')
> aluminum_pure.add_element('Al',1.0)
> aluminum_pure.set_density('g/cc', 2.7)
>
> aluminum_6061 = openmc.Material(1, 'aluminum 6061')
> aluminum_6061.add_element('Al', 0.9861)
> aluminum_6061.add_element('Mg', 0.008)
> aluminum_6061.add_element('Si', 0.004)
> aluminum_6061.add_element('Cu', 0.0015)
> aluminum_6061.add_element('Cr', 0.0004)
> aluminum_6061.set_density('g/cc', 2.7)
>
> germanium = openmc.Material(2, 'germanium')
> germanium.add_element('Ge', 1)
> germanium.set_density('g/cc', 5.35)
>
> iron = openmc.Material(3, 'iron')
> iron.add_element('Fe', 1)
> iron.set_density('g/cc', 7.874)
>
> air = openmc.Material(4, 'air')
> air.set_density('g/cc', 0.001205)
> air.add_element('N', 0.784431)
> air.add_element('O', 0.210748)
> air.add_element('Ar',0.0046)
>
> materials = openmc.Materials([aluminum_pure, aluminum_6061, germanium,iron,air])
> materials.export_to_xml()
>
> ### DEF SURFACES ###
>
> ge_radius_surf = openmc.ZCylinder(r = 7.0/2 - 0.1)
> endcap_radius_surf = openmc.ZCylinder(r = 7.0/2)
>
> ge_top_surf = openmc.ZPlane(z0=(13.4 - 1))
> endcap_top_surf = openmc.ZPlane(z0=(13.4))
>
> partition_surf = openmc.ZPlane(z0=0)
>
> bottom_flat_surf = openmc.ZPlane(z0=-11.6)
> bottom_rad_surf = openmc.ZCylinder(r=11.6/2)
>
> s = openmc.Sphere(r = 50, boundary_type='vacuum')
>
> ### DEF CELLS ###
>
> germanium_crystal = openmc.Cell(0,'crystal',germanium)
> germanium_crystal.region = -ge_radius_surf & -ge_top_surf & +partition_surf
>
> endcap = openmc.Cell(1, 'endcap', aluminum_pure)
> endcap.region = -endcap_radius_surf & -endcap_top_surf & +partition_surf & ~(-ge_radius_surf & -ge_top_surf)
>
> bottom_mount = openmc.Cell(2, 'mount', aluminum_6061)
> bottom_mount.region = -partition_surf & +bottom_flat_surf & -bottom_rad_surf
>
> surroundings = openmc.Cell(3, 'surroundings',fill=air)
> surroundings.region = -s
>
> ### DEF UNIVERSE ###
> universe = openmc.Universe(cells=[germanium_crystal,endcap,bottom_mount,surroundings])
> geometry = openmc.Geometry(universe)
> geometry.export_to_xml()
> universe.plot(width=(20.0, 50.0), origin=(0.0, 0.0, 5.0), basis='xz', colors=colors, color_by='material')
>
> ### DEF SOURCES ###
>
> source_1 = openmc.IndependentSource()
> source_1.space = openmc.stats.Point((0.0, 0.0, 30.0))
> source_1.angle = openmc.stats.Monodirectional([0.0, 0.0, -1.0])
> source_1.energy = openmc.stats.Discrete([800_000.0], [1.0])
> source_1.strength = 1e12
> source_1.particle = 'photon'
> source_1.domains = [germanium_crystal,endcap,bottom_mount,surroundings]
>
> openmc.config['chain_file'] = "./chain_simple.xml"
>
> source_Cs137 = openmc.IndependentSource()
> source_Cs137.space = openmc.stats.Point((0.0, 0.0, 30.0))
> source_Cs137.angle = openmc.stats.Isotropic()
> source_Cs137.energy = openmc.data.decay_photon_energy("Cs137")
> source_Cs137.particle = 'photon'
> source_Cs137.strength = np.sum(openmc.data.decay_photon_energy("Cs137").p) * 3e24 # Number of Atoms
> source_Cs137.domains = [germanium_crystal,endcap,bottom_mount,surroundings]
>
> ### DEF SETTINGS ###
>
> settings = openmc.Settings()
> settings.particles = 10**4
> settings.batches = 10
> settings.photon_transport = True
> settings.source = [source_1, source_Cs137]
> settings.run_mode = 'fixed source'
> settings.export_to_xml()
>
> ### DEF TALLIES ###
> tallies = openmc.Tallies()
>
> energy_bins = np.linspace(0, 1e6, 1001)
> energy_filter = openmc.EnergyFilter(energy_bins)
> cell_filter_ge = openmc.CellFilter(germanium_crystal)
> cell_filter_endcap = openmc.CellFilter(endcap)
>
> tally = openmc.Tally(name='pulse-height')
> tally.filters = [cell_filter_ge, energy_filter]
> tally.scores = ['pulse-height']
> tallies.append(tally)
> tallies.export_to_xml()
>
> ### RUN OPENMC ###
>
> openmc.run()
>
> ### POST PROCESSING ###
>
> sp = openmc.StatePoint('statepoint.10.h5')
> tally = sp.get_tally(name='pulse-height')
> pulse_height_values = tally.get_values(scores=['pulse-height']).flatten('C')
>
> energy_bin_centers = energy_bins[1:] + 0.5 * (energy_bins[1] - energy_bins[0])
>
> plt.figure()
> plt.semilogy(energy_bin_centers, pulse_height_values)
>
> # plot the strongest sources as vertical lines
> plt.axvline(x=800_000, color="red", ls=":") # source_1
> plt.axvline(x=661_700, color="red", ls=":") # source_2
>
> plt.xlabel('Energy [eV]')
> plt.ylabel('Counts')
> plt.title('Pulse Height Values')
> plt.grid(True)
> plt.tight_layout()
>
> ### CHECK PULSE HEIGHT VARIABLE ###
>
> print(pulse_height_values)
>
> sp.close()
