Energy Spectrum Attenuation

[MadGander]

TL;DR

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...

 

[Mark44]

Should be relatively simple to debug. I count only 9 lines or so of executable statements with all the rest being data.

 

[Vanadium 50]

all the rest being data.

And a typo in data can never cause a bug? You've lived a charmed life, my friend!

I'm not even sure there is a problem. If you filter out the low energy stuff, the mean energy moves up, not down. Maybe its OK. Maybe irs not.

 

[Mark44]

And a typo in data can never cause a bug?

Of course it can, but that's not something that we as observers of the code can find. What I meant by my comment was that with so few lines of executable code, the problem isn't there. By implication, it's very likely a problem with the data.

 

[berkeman]

I guess you could move it over to the Nuke section.

Update -- Moved to the Nuclear Engineering forum. @MadGander -- can you address some of the questions raised so far for you? Thanks.

 

[Astronuc]

TL;DR Summary: Higher average energy after passing through lead

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, . . .

Yes, something is wrong if one calculates an 'attenuated' energy of 74 keV from an initial X-ray (photon) energy of 33 keV. One may see a dose build up inside the surface as electrons and photons interact with the solid, and the photons/electrons lose energy as they scatter on the atomic electrons. The dose is just the integration of all the particles (photons and electrons) losing energy.

One should be careful in how one sums (integrates) the energy in a unit volume. A photon can lose no energy (Rayleigh scattering, or coherent or elastic scattering), some energy (Compton scattering, or inelastic scattering), or all of the energy (photoelectric effect).

https://tech.snmjournals.org/content/33/1/3

At some distance into the solid, some fraction (x) will have interacted with electrons (scattering) while the other fraction (1-x) have not yet interacted/scatter. Also, a lot depends on which electron in an atom (for Pb, Z = 82) the electron interacts.

Is one applying the physics correctly?

 

[Alex A]

The context is in @MadGander's other thread. Consider the aluminium filter on an X-ray machine. Those low energy X-rays from the bare tube don't contribute to the image because they don't make it through the soft tissue, causing sunburn like damage. The filter makes the beam 'harder' - the average energy of a photon making it through the filter is increased. No photons have gained energy going through the filter. It's just that most low energy photons are stopped dead, and higher energy photons are barely attenuated. What must reduce to satisfy conservation of energy is the total power in the beam.

 

[mfb]

It's just that most low energy photons are stopped dead, and higher energy photons are barely attenuated.

This is the right answer. The shielding doesn't change any photon energy, it just changes the distribution by removing some of the photons. More shielding would raise the average photon energy even more while further reducing the intensity.

 

[Vanadium 50]

This possibility was brought up in Post #9. No response from the OP.