# coupled photon transport can be enabled for a neutron transport
# with the settings.photon_transport = True
# that is fairly standard and we need to enable photon production and transport
# for tallies like heat, dose and it is important to transport the photon
# as it will not deposit its energy locally.
# Charged particles will deposit their energy locally so they don't need to be transported

# However this example is somewhat different
# This examples shows that a 4MeV photon source creates neutrons
# We do this by creating a photon source and tallying neutron flux

# You could try changing the material to a lower Z atom like Be or Li and the
# flux should go up partly because these shield photons to a lesser extent.

import openmc

# MATERIALS
material = openmc.Material()
material.add_element('Fe', 1)
material.set_density('g/cm3', 7)
my_materials = openmc.Materials([material])

# GEOMETRY
# surfaces
vessel_inner = openmc.Sphere(r=500, boundary_type='vacuum')
# cells
inner_vessel_region = -vessel_inner
inner_vessel_cell = openmc.Cell(region=inner_vessel_region)
inner_vessel_cell.fill = material
my_geometry = openmc.Geometry([inner_vessel_cell])

# SIMULATION SETTINGS
# Instantiate a Settings object
my_settings = openmc.Settings()
my_settings.batches = 10
my_settings.inactive = 0
my_settings.particles = 500
my_settings.run_mode = 'fixed source'
# Create a photon point source
my_source = openmc.Source()
my_source.space = openmc.stats.Point((0, 0, 0))
my_source.angle = openmc.stats.Isotropic()
my_source.particle = 'photon'  # default would otherwise be neutron
my_source.energy = openmc.stats.Discrete([4e6], [1])
my_settings.source = my_source
# added a cell tally for neutron flux production
cell_filter = openmc.CellFilter(inner_vessel_cell)
flux_tally = openmc.Tally(name='flux')
flux_tally.filters = [cell_filter]
flux_tally.scores = ['flux']
my_tallies = openmc.Tallies([flux_tally])

# Run OpenMC!
model = openmc.model.Model(my_geometry, my_materials, my_settings, my_tallies)
sp_filename = model.run()

# open the results file
sp = openmc.StatePoint(sp_filename)

# access the tally using pandas dataframes
flux_tally = sp.get_tally(name='flux')

# printing the tally shows that we get a neutron flux as the tally is not 0
print(flux_tally.mean, flux_tally.std_dev)
# printed results are ...
# >>> [[[14.42136056]]] [[[0.14143424]]]