Energy Spectrum Attenuation

[MadGander]

TL;DR Summary

Higher average energy after passing through lead

Hello,

I'm currently trying to compare theoretical results with an MCNP simulation. I'm using two discrete sets of data, intensity (probability) and linear attenuation coefficient, both functions of energy, to produce an attenuated energy spectrum after x-rays have passed through a thin layer of lead. I've been running through the calculations and I'm getting a higher average attenuated energy (~74 keV) than initial average energy (~33 keV). My guess is I'm doing something wrong somewhere in my calculations, but I've checked and re-checked to the upteenth time and can't find any errors. I'll attach my python code below. Appreciate any assistance because I'm literally going crazy.

Thanks

 

[MadGander]

Here is the Python code

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
import math


# Lead

Energies = np.linspace(1,100,199)
LAC = np.array([
    5.87E+04, 2.65E+04, 1.44E+04, 1.65E+04, 2.21E+04,
    1.74E+04, 1.40E+04, 1.05E+04, 8.16E+03, 6.45E+03,
    5.20E+03, 4.25E+03, 3.52E+03, 2.96E+03, 2.52E+03,
    2.15E+03, 1.86E+03, 1.62E+03, 1.42E+03, 1.25E+03,
    1.11E+03, 9.88E+02, 8.85E+02, 7.96E+02, 7.19E+02,
    1.63E+03, 1.47E+03, 1.34E+03, 1.22E+03, 1.57E+03,
    1.68E+03, 1.55E+03, 1.44E+03, 1.33E+03, 1.24E+03,
    1.16E+03, 1.08E+03, 1.01E+03, 9.50E+02, 8.90E+02,
    8.35E+02, 7.85E+02, 7.39E+02, 6.97E+02, 6.57E+02,
    6.21E+02, 5.88E+02, 5.57E+02, 5.28E+02, 5.01E+02,
    4.76E+02, 4.53E+02, 4.31E+02, 4.11E+02, 3.92E+02,
    3.74E+02, 3.57E+02, 3.42E+02, 3.27E+02, 3.13E+02,
    3.00E+02, 2.87E+02, 2.75E+02, 2.64E+02, 2.54E+02,
    2.44E+02, 2.34E+02, 2.25E+02, 2.17E+02, 2.09E+02,
    2.01E+02, 1.94E+02, 1.87E+02, 1.80E+02, 1.74E+02,
    1.68E+02, 1.62E+02, 1.57E+02, 1.52E+02, 1.47E+02,
    1.42E+02, 1.38E+02, 1.33E+02, 1.29E+02, 1.25E+02,
    1.21E+02, 1.18E+02, 1.14E+02, 1.11E+02, 1.07E+02,
    1.04E+02, 1.01E+02, 9.85E+01, 9.57E+01, 9.31E+01,
    9.06E+01, 8.81E+01, 8.57E+01, 8.35E+01, 8.13E+01,
    7.91E+01, 7.71E+01, 7.51E+01, 7.32E+01, 7.14E+01,
    6.96E+01, 6.79E+01, 6.62E+01, 6.46E+01, 6.31E+01,
    6.16E+01, 6.01E+01, 5.87E+01, 5.73E+01, 5.60E+01,
    5.48E+01, 5.35E+01, 5.23E+01, 5.12E+01, 5.01E+01,
    4.90E+01, 4.79E+01, 4.69E+01, 4.59E+01, 4.49E+01,
    4.40E+01, 4.31E+01, 4.22E+01, 4.13E+01, 4.05E+01,
    3.97E+01, 3.89E+01, 3.81E+01, 3.74E+01, 3.66E+01,
    3.59E+01, 3.52E+01, 3.46E+01, 3.39E+01, 3.33E+01,
    3.27E+01, 3.21E+01, 3.15E+01, 3.09E+01, 3.04E+01,
    2.98E+01, 2.93E+01, 2.88E+01, 2.83E+01, 2.78E+01,
    2.73E+01, 2.68E+01, 2.64E+01, 2.59E+01, 2.55E+01,
    2.51E+01, 2.47E+01, 2.43E+01, 2.39E+01, 2.35E+01,
    2.31E+01, 2.27E+01, 2.24E+01, 2.20E+01, 2.17E+01,
    2.13E+01, 2.10E+01, 2.07E+01, 2.04E+01, 2.01E+01,
    1.98E+01, 1.95E+01, 1.92E+01, 1.89E+01, 1.86E+01,
    8.26E+01, 8.15E+01, 8.03E+01, 7.91E+01, 7.80E+01,
    7.69E+01, 7.58E+01, 7.47E+01, 7.37E+01, 7.27E+01,
    7.17E+01, 7.07E+01, 6.98E+01, 6.88E+01, 6.79E+01,
    6.70E+01, 6.61E+01, 6.52E+01, 6.44E+01, 6.35E+01,
    6.27E+01, 6.19E+01, 6.11E+01, 6.03E+01
])

I_0 = np.linspace(1,100,199)
P_0 = np.array([   
    0,
    0,
    5.11347E-06,
    0.000121555,
    3.08604E-05,
    0,
    0,
    8.31027E-06,
    8.7771E-05,
    0.000359376,
    0.000946079,
    0.001877169,
    0.003548272,
    0.005150231,
    0.006742632,
    0.016790123,
    0.014374974,
    0.012801486,
    0.024634899,
    0.011637658,
    0.015204991,
    0.0107191,
    0.008945965,
    0.009214431,
    0.015634477,
    0.009179067,
    0.009152242,
    0.00943935,
    0.010478518,
    0.010235272,
    0.010345676,
    0.010531003,
    0.010234591,
    0.01073557,
    0.010859254,
    0.011212683,
    0.011140717,
    0.011089523,
    0.010867856,
    0.010987345,
    0.011155824,
    0.010809004,
    0.011037909,
    0.010752565,
    0.010937409,
    0.010729171,
    0.010745011,
    0.010583666,
    0.010481057,
    0.010277613,
    0.010390617,
    0.010158208,
    0.010171552,
    0.009960702,
    0.010035846,
    0.009819078,
    0.009828614,
    0.009517127,
    0.009333436,
    0.009406346,
    0.009269727,
    0.008867548,
    0.008940091,
    0.008521914,
    0.008432146,
    0.008507657,
    0.008371625,
    0.008354725,
    0.00794921,
    0.007963582,
    0.007768016,
    0.007666614,
    0.00774671,
    0.007514888,
    0.007287714,
    0.007189879,
    0.007027001,
    0.006843111,
    0.006814986,
    0.006970499,
    0.006564386,
    0.006477514,
    0.006288494,
    0.006329008,
    0.006155955,
    0.006046507,
    0.005918102,
    0.005944737,
    0.00578743,
    0.005658994,
    0.00563903,
    0.005555042,
    0.005472806,
    0.005113881,
    0.005375275,
    0.005027333,
    0.005019937,
    0.005089112,
    0.004891952,
    0.004766316,
    0.004586465,
    0.004527969,
    0.004351737,
    0.004400257,
    0.004311621,
    0.004322532,
    0.004149311,
    0.004076936,
    0.003974694,
    0.004101935,
    0.003906201,
    0.003932123,
    0.003898879,
    0.003951311,
    0.003648835,
    0.008549116,
    0.003395487,
    0.012448897,
    0.003362305,
    0.003395445,
    0.003248954,
    0.003229609,
    0.003147635,
    0.002916915,
    0.002987517,
    0.002910967,
    0.002877974,
    0.00273077,
    0.002851213,
    0.002961867,
    0.002615342,
    0.002670858,
    0.002596836,
    0.005511989,
    0.002374152,
    0.002347264,
    0.002388104,
    0.002970019,
    0.002251517,
    0.001893168,
    0.001914275,
    0.002046383,
    0.001923139,
    0.001900081,
    0.00176021,
    0.002097179,
    0.001679893,
    0.001536802,
    0.001690877,
    0.002112264,
    0.001642358,
    0.001356783,
    0.00140141,
    0.001483782,
    0.001486027,
    0.001390132,
    0.001416076,
    0.001205739,
    0.001301917,
    0.001234389,
    0.001154041,
    0.001219586,
    0.00121325,
    0.001145995,
    0.001201994,
    0.001040079,
    0.001012389,
    0.000876481,
    0.001004057,
    0.000925759,
    0.000918481,
    0.000809436,
    0.000821415,
    0.000789026,
    0.000813591,
    0.00068724,
    0.000703664,
    0.000756499,
    0.000682395,
    0.000586809,
    0.000578912,
    0.00049449,
    0.000509462,
    0.000464834,
    0.00045288,
    0.000404728,
    0.000394017,
    0.000356244,
    0.000352201,
    0.000281882,
    0.000271873,
    0.00022954,
    0.000237268,
    0.000188778,
    0.000170614,
    0.000129775,
    9.31903E-05,
    6.61896E-05,
    2.2146E-05
])

thickness = 0.08

P_final = np.zeros(199)

for i in range(len(P_0)):
   P_final[i] = P_0[i] * np.exp(-LAC[i] * thickness)
    
print(P_final)



# Calculate mean energies based on weighted averages
MeanFinalE = np.sum(Energies * P_final) / np.sum(P_final) if np.sum(P_final) > 0 else 0
MeanInitialE = np.sum(Energies * P_0) / np.sum(P_0) if np.sum(P_0) > 0 else 0

# Output results
print(f"Mean Initial Energy: {MeanInitialE:.2f} keV")
print(f"Mean Final Energy: {MeanFinalE:.2f} keV")

 

[berkeman]

Welcome to PF.

Is this mostly a Python question, or an MCNP question? MCNP questions are usually posted in our Nuclear Engineering forum, where folks have expertise in that software. Let me know if you'd like the thread moved, or if it is okay to stay in this forum.

 

[MadGander]

Welcome to PF.

Is this mostly a Python question, or an MCNP question? MCNP questions are usually posted in our Nuclear Engineering forum, where folks have expertise in that software. Let me know if you'd like the thread moved, or if it is okay to stay in this forum.

My results were calculated using Python and don't directly pertain to MCNP. Wasn't sure exactly where to post this question, but I guess you could move it over to the Nuke section. Thanks.

 

[berkeman]

We can leave it here for now and see how it goes.

 

[Vanadium 50]

"Debug my code for me" is a pretty heavy lift. Especially when a) you don't show what you have already done to debug it, b) have pared it down to the bare minimum to show the problems, and c) the code is almost completely uncommented. "Here are hundreds of numbers - you guys tell me which one is wrong" is, as I said, a big ask.

 

[berkeman]

"Debug my code for me" is a pretty heavy lift. Especially when a) you don't show what you have already done to debug it, b) have pared it down to the bare minimum to show the problems, and c) the code is almost completely uncommented. "Here are hundreds of numbers - you guys tell me which one is wrong" is, as I said, a big ask.

All valid points.

@MadGander -- please help us out here...