Cs3.5MeV spectrum

User Support
1

Hi all,
I am trying to get the escape peaks for a 3.5MeV Cs source. I don’t know where I am doing the mistake.
The input is given below:

import openmc
import numpy as np
import matplotlib.pyplot as plt
import math
import openmc_source_plotter

# MATERIALS

# Tungsten is a very good photon shield, partly due to its high Z number and electrons
detector_material = openmc.Material()
detector_material.add_element("Ge",1)
#detector_material.add_nuclide("Ge70", 0.2052, percent_type="ao")
#detector_material.add_nuclide("Ge72", 0.2745, percent_type="ao")
#detector_material.add_nuclide("Ge73", 0.0776, percent_type="ao")
#detector_material.add_nuclide("Ge74", 0.3670, percent_type="ao")
#detector_material.add_nuclide("Ge76", 0.0775, percent_type="ao")
detector_material.set_density('g/cm3', 5.323)

# Air
air = openmc.Material(name='Air')
air.set_density('g/cm3', 0.0012050)
air.add_element('N', 0.755268, 'wo')
air.add_element('O', 0.231781, 'wo')
air.add_element('Ar', 0.012827, 'wo')
air.add_nuclide('C0', 0.000124, 'wo')
air.remove_nuclide('O18')

material = openmc.Materials([detector_material,air])

sphere_surface = openmc.Sphere(r=50, boundary_type="vacuum")

detector_outer_surface=openmc.Sphere(r=15)
detector_inner_surface=openmc.Sphere(r=1)
detector_region=(-detector_outer_surface & +detector_inner_surface)
#detector_cylinder=openmc.ZCylinder(r=8)
#detector_region=(-detector_cylinder & -sphere_surface)

detector_cell = openmc.Cell(region=detector_region)
detector_cell.fill = detector_material

void_space_cell = openmc.Cell(region=-sphere_surface & ~detector_region)
void_space_cell.fill= air

universe = openmc.Universe(cells=[detector_cell,void_space_cell])
geometry = openmc.Geometry(universe)

source = openmc.Source()
source.space = openmc.stats.Point((0, 0, 20))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([3.5e6], [1])
#source.energy = openmc.stats.Discrete([1.1732e6,1.3325e6], [0.5, 0.5])
#source.energy = openmc.stats.Discrete([0.05316e6, 0.07961e6 , 0.08099e6 ,0.2763e6 ,0.30285e6 ,0.35601e6 ,0.38384e6], [0.0214, 0.0265 ,0.3209 ,0.0716 ,0.1834 ,0.6205 ,0.0894])
source.particle='photon'

setting = openmc.Settings()
setting.batches = 2
setting.inactive = 0 # the default is 10, which would be wasted computing for us
setting.particles = 1000000
setting.run_mode = 'fixed source'
setting.photon_transport = True  # This line is required to switch on photons tracking
setting.source = source

tallies = openmc.Tallies()

# sets up filters for the tallies
particle_filter = openmc.ParticleFilter(['photon'])  # note the use of photons here
#energy_filter = openmc.EnergyFilter(np.linspace(0,4e6,1958))


# setup the filters for the cell tally
#cell_filter = openmc.CellFilter(detector_cell,void_space_cell) 

# create the tally
cell_spectra_tally = openmc.Tally(name='cell_spectra_tally')
cell_spectra_tally.scores = ['flux']
cell_spectra_tally.filters = [particle_filter,cell_filter] #energy_filter
tallies.append(cell_spectra_tally)

# combine all the required parts to make a model
model = openmc.model.Model(geometry, material, setting, tallies)

# remove old files and runs OpenMC
!rm *.h5
results_filename = model.run()

results = openmc.StatePoint(results_filename)

#extracts the tally values from the simulation results
cell_tally = results.get_tally(name='cell_spectra_tally')

# imports a convenient plotting package
from spectrum_plotter import plot_spectrum_from_tally

plot_spectrum_from_tally(
    spectrum={'cell tally':cell_tally},
    x_label = "energy [eV]",
    y_label = "Photons per cm2",
    plotting_package="plotly",  # matplotlib is another option, but plotly option provides a dropdown menu for changing axis scale
)
2

@paulromano Can you please provide your comments

3

If you’re running into problems with openmc_source_plotter, that package is not officially maintained by the OpenMC development team. Perhaps @Shimwell can comment on this, or you can create an issue on the following repository:

4

@asadullahamin I am not sure what error message you were getting can you please let us know the error message otherwise it makes it hard to reproduce the same bug.

However I did have a go at running your script. I had to change

cell_filter = openmc.CellFilter(detector_cell,void_space_cell)

to

cell_filter = openmc.CellFilter([detector_cell,void_space_cell])

to get it to run and produce an error.

Was this you error?

 raise ValueError("EnergyFilter was not found in spectra tally")
ValueError: EnergyFilter was not found in spectra tally

If so then this indicates that you need an EnergyFilter in the tally to get a spectrum plot. This makes sense to me as you need to split up the tally into energy bins to plot a graph with energy on the axis

Here is an example of a tally with an energy filter

However please do raise an issue on the repo that Paul linked if you have more details

5

Hi all, I am getting the spectrum plots but I am not getting the escape peaks. I am not sure whether pulse height spectra can be plotted with spectrum plots.
I have used the energy filter.
Do I have to use pulse height tally to get that plots?
I am not sure whether it is available with the current version?

6

I have tried using energy filter, you can see it there but it is commented. I tried running with that to.

7

Pulse height tallies are being proposed on this pull request but are not quite in OpenMC yet