Hello,
I have been trying to simulate a 2"x2" NaI detector and Cs137 point source. The goal is to obtain the same amount of counts under the curve as the real measurements. However, my simulation results are too high (no consistent scaling factor).
For the real case, I used Genie 2K and obtained the counts by setting a region of interest over the entire curve. The background was 46 cpm and dead time was negligible (<1%).
The source used was a Cs-137 check source of initially 9383 Bq (July 1 2013), measured to be 7215 Bq at the time of the experiment.
I then placed the source at distances of 2.5, 5, and 10 cm from the detector, which was set standing up on the table. A 0 cm measurement was taken with the check source placed underneath the detector. Note that this value I am less concerned over, as the point source in the simulation has been placed directly against the detector, whereas the real check source has some thickness. Thus, the results for this case may not be representative of the real case.
V0.15.0 of OpenMC was used. The OpenMC model (pasted below) is kept as simple as can be, modelling only air, the detector crystal, and a thin aluminum casing on the walls only (with a large, r=20, out-of-bounds sphere surrounding it all). The source is an isotropic point source moved along the X-axis away from the detector wall.
The source strength is found in the manner used in the gamma detector example, with the mass set to ~2.24 ng (found by comparing the specific activity of 3.2 TBq/g to the known activity of 7215 Bq). At runtime, this calculates a strength of 6665 Bq.
The tally uses an energy filter as described in the same example, with the 0-1 keV bin cut out. The “pulse-height” score is used and the counts obtained from the get_reshaped_data sum. I am not sure if this is the correct method, as in this portion I’m aping what I’ve seen in other examples.
The results are as follows:
| Case | Net Real (cpm) | Simulated (cpm) |
|---|---|---|
| Contact | 239 | 659 |
| 2.5 cm | 72 | 92 |
| 5 cm | 24 | 46 |
| 10 cm | 12 | 18 |
The main questions I have are:
- Am I misusing the pulse-height tally? If so, is there a more appropriate tally for comparing the counts?
- Am I retrieving the tally data correctly for purposes of finding counts?
- Are there any limits of OpenMC that I need to consider in post-processing, e.g., intrinsic efficiency (as there can be no modelling of detector electronics)?
Thank you for your time.
CODE
# Detector: 2x2 NaI
# Source: Cs-137, 7215 Bq (originally 9383 Bq on July 1 2013)
# OpenMC Version 0.15.0
# IMPORTS
import openmc as omc
import numpy as np
from scipy.constants import Avogadro
# FUNCTIONS
def write_strings_to_txt(stringlist, filename, dir_dest, append=False):
mode = 'w'
if append:
mode = 'a'
with open(dir_dest + filename, mode) as f:
f.write('\n'.join(stringlist))
# MATERIALS
FRACTION_ATOMIC = "ao"
mat_air = omc.Material(name="air_dry_near_sea_level_noC")
mat_air.set_density("g/cm3", 0.001205)
mat_air.add_element('C', 0.000150, percent_type=FRACTION_ATOMIC)
mat_air.add_element('N', 0.784431, percent_type=FRACTION_ATOMIC)
mat_air.add_element('O', 0.210748, percent_type=FRACTION_ATOMIC)
mat_air.add_element('Ar', 0.004671, percent_type=FRACTION_ATOMIC)
mat_NaI = omc.Material(name="NaI")
mat_NaI.set_density("g/cm3", 3.667)
mat_NaI.add_element('Na', 0.499999, percent_type=FRACTION_ATOMIC)
mat_NaI.add_element('I', 0.500001, percent_type=FRACTION_ATOMIC)
aluminum = omc.Material(name="aluminum")
aluminum.add_element('Al', 1)
aluminum.set_density('g/cm3', 2.7)
materials = omc.Materials([mat_air, mat_NaI, aluminum])
materials.cross_sections = r"/home/USER/endfb-viii.0-hdf5/cross_sections.xml"
materials.export_to_xml()
# MEASUREMENTS
IN_TO_CM = 2.54
meas_detector_radius = 1 * IN_TO_CM
meas_detector_height = 2 * IN_TO_CM
# POSITIONS
posz_detector_bottom = 0
posz_detector_top = posz_detector_bottom + meas_detector_height
# GEOMETRY
surf_detector_cyl = omc.ZCylinder(r=meas_detector_radius, x0=0, y0=0)
surf_detector_bottom = omc.ZPlane(z0=posz_detector_bottom)
surf_detector_top = omc.ZPlane(z0=posz_detector_top)
surf_oob_sphere = omc.Sphere(r=20, x0=0, y0=0, z0=0, boundary_type="vacuum")
surf_casing_cyl_outer = omc.ZCylinder(r=meas_detector_radius+0.05, x0=0, y0=0)
region_detector = -surf_detector_cyl & +surf_detector_bottom & -surf_detector_top
region_casing = +surf_detector_cyl & -surf_casing_cyl_outer & +surf_detector_bottom & -surf_detector_top
region_air = -surf_oob_sphere & ~region_detector & ~region_casing
region_oob = +surf_oob_sphere
cell_detector = omc.Cell(fill=mat_NaI, region=region_detector, name="TARGET")
cell_casing = omc.Cell(fill=aluminum, region=region_casing, name="casing")
cell_air = omc.Cell(fill=mat_air, region=region_air)
cell_oob = omc.Cell(region=region_oob)
cell_list = [cell_detector, cell_casing, cell_air, cell_oob]
universe = omc.Universe(cells=cell_list)
geometry = omc.Geometry(universe)
geometry.export_to_xml()
# Plot Geometry
plot_width = 20
plot_pixel_resolution = (400, 400)
view_xy = "xy"
view_xz = "xz"
view_yz = "yz"
viewlist = [view_xy, view_xz, view_yz]
plot_color_by_material = "material"
DictColors = {mat_air: (255, 255, 255), mat_NaI: (177, 156, 217), aluminum: (100, 100, 100)}
plotslist = []
for viewindx in range(len(viewlist)):
plotname = "geometryview_" + viewlist[viewindx]
plot = omc.Plot(name=plotname)
plot.width = (plot_width, plot_width)
plot.pixels = plot_pixel_resolution
plot.color_by = plot_color_by_material
plot.colors = DictColors
plot.basis = viewlist[viewindx]
plot.filename = plotname
plot.origin = (0, 0, 0)
plotslist.append(plot)
# Export Plots
plots = omc.Plots(plotslist)
plots.export_to_xml()
# Create Images
omc.plot_geometry()
# SIMULATE
distances = [0, 2.5+meas_detector_radius, 5+meas_detector_radius, 10+meas_detector_radius] # Placed outside surface
for x in range(len(distances)):
meas_source_distance = distances[x]
omc.config['chain_file'] = "/home/USER/PycharmProjects/Cs137_PointSource_Test/chain_simple.xml"
energy_bins = np.linspace(0, 1e6, 1001)[1:] # Cut off the first bin to avoid no-energy cases
half_life_Cs137_seconds = 30.17 * 365 * 24 * 60 * 60
half_life_Ba137_m1_seconds = 2.55 * 60
branching_ratio = 0.946
fraction_atoms_Ba137_m1 = (half_life_Ba137_m1_seconds * branching_ratio) / half_life_Cs137_seconds
number_atoms_Ba137_m1 = (2.244168*10**-9 / 137) * Avogadro * fraction_atoms_Ba137_m1
number_atoms_Cs137 = (2.244168*10**-9 / 137) * Avogadro * (1 - fraction_atoms_Ba137_m1)
cs137_strength = np.sum(omc.data.decay_photon_energy("Cs137").p) * number_atoms_Cs137
ba137_strength = np.sum(omc.data.decay_photon_energy("Ba137_m1").p) * number_atoms_Ba137_m1 # 7264.53 > Secular equilibrium
sourceBa = omc.IndependentSource(angle=omc.stats.Isotropic())
sourceBa.space = omc.stats.Point((meas_source_distance, 0, 0))
sourceBa.energy = omc.data.decay_photon_energy("Ba137_m1")
sourceBa.strength = ba137_strength
sourceBa.particle = "photon"
sourceCs = omc.IndependentSource(angle=omc.stats.Isotropic())
sourceCs.space = omc.stats.Point((meas_source_distance, 0, 0))
sourceCs.energy = omc.data.decay_photon_energy("Cs137")
sourceCs.strength = cs137_strength
sourceCs.particle = "photon"
settings = omc.Settings()
settings.batches = 5
settings.particles = 10**5
settings.inactive = 0
settings.photon_transport = True
settings.run_mode = "fixed source"
settings.source = [sourceBa, sourceCs]
cell_filter = omc.CellFilter(cell_detector)
energy_filter = omc.EnergyFilter(energy_bins)
# photon_particle_filter = omc.ParticleFilter("photon")
name_tally = "pulse-height"
cell_tally = omc.Tally(name=name_tally)
cell_tally.filters = [cell_filter, energy_filter]
cell_tally.scores = [name_tally]
tallies = omc.Tallies([cell_tally])
tallies.export_to_xml()
model = omc.Model(geometry, materials, settings, tallies)
statepoint_filename = model.run()
S_TO_MIN = 60
with omc.StatePoint(statepoint_filename) as statepoint:
tally_result = statepoint.get_tally(name=name_tally)
sum_tally = tally_result.get_reshaped_data(value="sum").flatten()[0] * S_TO_MIN
std_dev = tally_result.get_reshaped_data(value="std_dev").flatten()[0] * S_TO_MIN
rel_err = tally_result.get_reshaped_data(value="rel_err").flatten()[0] * S_TO_MIN
output = ["Score Type: " + name_tally,
f"Source Distance of (cm, from origin): {meas_source_distance:.1f}",
f"Tally (cpm): {sum_tally:.1f} +/- {rel_err:.1f}",
'\n']
path = r"/home/USER/PycharmProjects/Cs137_PointSource_Test/"
name = "output" + name_tally + ".txt"
write_strings_to_txt(output, name, path, append=True)
write_strings_to_txt(["----------------------------------\n"], name, path, append=True)