tail_tau needs to be unreasonable large for lines in MeV range

[joefowler]

Original report by Galen O'Neil (Bitbucket: oneilg, GitHub: oneilg).

---

For True Bq we’re trying to fit the spectrum with see with tails. The tails, by eye, are on the 10-100 keV length scale. But when I use a tail_tau parameter of that size the spectrum is unchanged by eye. Instead I have to use a tail on the length scale of 1e4 keV to get a reasonable fit. Seems like the length scale for tails is energy dependent? It shouldn’t be.

‌

Here is an example of the spectrum we want to fit, with my most successful fit so far.

‌

‌

Here is the output of the below minimum reproducing example:

‌

‌

# -*- coding: utf-8 -*-

import mass
import numpy as np
import pylab as plt
plt.ion()

fwhm = 5000.0
tail_tau_big = 1e7
tail_tau_reasonable = 5e4
dph_de=1
def cleanspectrum_fn(x): return mass.getmodel(5637.82e3).spect.pdf(x, instrument_gaussian_fwhm=fwhm)
bin_width=1e3
energy=np.arange(5.5e6, 5.8e6, bin_width)
length_lo_big = tail_tau_big*dph_de/bin_width
length_lo_reasonable = tail_tau_reasonable*dph_de/bin_width
length_hi = 0
peak_ph = 5637.82e3
integral = 47466.4
background = 28.6
tail_frac = 0.5
tail_frac_hi = 0
spect_bigtau = mass.line_fits._smear_exponential_tail(
cleanspectrum_fn, energy, fwhm, tail_frac, length_lo_big, tail_frac_hi, length_hi)
spect_reasonabletau = mass.line_fits._smear_exponential_tail(
cleanspectrum_fn, energy, fwhm, tail_frac, length_lo_reasonable, tail_frac_hi, length_hi)
plt.figure()
plt.plot(energy, spect_bigtau, label=f"tail_tau={tail_tau_big:.2e}")
plt.plot(energy, spect_reasonabletau, label=f"tail_tau={tail_tau_reasonable:.2e}")
plt.legend()
plt.xlabel("energy (eV)")
plt.ylabel("counts  (arb)")

‌

[joefowler]

Original comment by Joseph Fowler (Bitbucket: joe_fowler, ).

---

Weird. That doesn’t seem right. I can’t think of what the problem would be??

[joefowler]

Original comment by Galen O'Neil (Bitbucket: oneilg, GitHub: oneilg).

---

I think I must not be doing the MRE right, because I added a scale factor to bring it down to “normal x-ray energy range” and the results are the same. The reasonable tail is not noticeable, while the extra large tail tau seems about normal.

‌

# -*- coding: utf-8 -*-

import mass
import numpy as np
import pylab as plt
plt.ion()

scale = 3e-4
fwhm = 5000.0*scale
tail_tau_big = 1e7*scale
tail_tau_reasonable = 5e4*scale
dph_de=1
def cleanspectrum_fn(x): return mass.getmodel(5637.82e3*scale).spect.pdf(x, instrument_gaussian_fwhm=fwhm)
bin_width=1e3
energy=np.arange(5.5e6, 5.8e6, bin_width)*scale
length_lo_big = tail_tau_big*dph_de/bin_width
length_lo_reasonable = tail_tau_reasonable*dph_de/bin_width
length_hi = 0
peak_ph = 5637.82e3*scale
integral = 47466.4
background = 28.6
tail_frac = 0.5
tail_frac_hi = 0
spect_bigtau = mass.line_fits._smear_exponential_tail(
cleanspectrum_fn, energy, fwhm, tail_frac, length_lo_big, tail_frac_hi, length_hi)
spect_reasonabletau = mass.line_fits._smear_exponential_tail(
cleanspectrum_fn, energy, fwhm, tail_frac, length_lo_reasonable, tail_frac_hi, length_hi)
plt.figure()
plt.plot(energy, spect_bigtau, label=f"tail_tau={tail_tau_big:.2e}")
plt.plot(energy, spect_reasonabletau, label=f"tail_tau={tail_tau_reasonable}")
plt.legend()
plt.xlabel("energy (eV)")
plt.ylabel("counts  (arb)")

‌

[joefowler]

Original comment by Galen O'Neil (Bitbucket: oneilg, GitHub: oneilg).

---

Ok I made a better minimum reproducing example to show the problem. It makes 3 plots

1. MnKalpha lineshape with and without tail (looks as expected)
2. A monochromatic line at the MnKAlpha energy with and without a tail (looks as expected)
3. Same as 2, but every energy involved is multiplied by scale=1e3 (no tail is visible!).

‌

import numpy as np
import mass
import pylab as plt
plt.ion()

line = mass.spectra["MnKAlpha"]
model = line.model(has_tails=True)
params = model.make_params()
params["fwhm"].set(2)
params["tail_tau"].set(30)
params["peak_ph"].set(mass.STANDARD_FEATURES["MnKAlpha"])
params["tail_frac"].set(0.5)
bin_centers = np.arange(5800,6000,1)
counts_tail = model.eval(params=params, bin_centers=bin_centers)

params["tail_frac"].set(0.0)
counts_notail = model.eval(params=params, bin_centers=bin_centers)

plt.figure()
plt.plot(bin_centers, counts_tail, label="counts_tail")
plt.plot(bin_centers, counts_notail, label="counts_notail")
plt.xlabel("energy (eV)")
plt.ylabel("counts (arb)")
plt.title("MnKAlpha")
plt.legend()

_model = mass.getmodel(mass.STANDARD_FEATURES["MnKAlpha"])
model=_model.spect.model(has_tails=True)
params = model.make_params()
params["fwhm"].set(2)
params["tail_tau"].set(30)
params["peak_ph"].set(mass.STANDARD_FEATURES["MnKAlpha"])
params["tail_frac"].set(0.5)
bin_centers = np.arange(5800,6000,1)
counts_tail = model.eval(params=params, bin_centers=bin_centers)

params["tail_frac"].set(0.0)
counts_notail = model.eval(params=params, bin_centers=bin_centers)

plt.figure()
plt.plot(bin_centers, counts_tail, label="counts_tail")
plt.plot(bin_centers, counts_notail, label="counts_notail")
plt.xlabel("energy (eV)")
plt.ylabel("counts (arb)")
plt.title("MnKAlpha1 mono")
plt.legend()

scale = 1e3
_model = mass.getmodel(mass.STANDARD_FEATURES["MnKAlpha"]*scale)
model=_model.spect.model(has_tails=True)
params = model.make_params()
params["fwhm"].set(2*scale)
params["tail_tau"].set(30*scale)
params["peak_ph"].set(mass.STANDARD_FEATURES["MnKAlpha"]*scale)
params["tail_frac"].set(0.5)
bin_centers = np.arange(5800,6000,1)*scale
counts_tail = model.eval(params=params, bin_centers=bin_centers)

params["tail_frac"].set(0.0)
counts_notail = model.eval(params=params, bin_centers=bin_centers)

plt.figure()
plt.plot(bin_centers, counts_tail, label="counts_tail")
plt.plot(bin_centers, counts_notail, label="counts_notail")
plt.xlabel("energy (eV)")
plt.ylabel("counts (arb)")
plt.title(f"MnKAlpha1 mono * {scale=}")
plt.legend()

‌

[joefowler]

Original comment by Galen O'Neil (Bitbucket: oneilg, GitHub: oneilg).

---

Ok I made another MRE that works just with _smear_exponential_tail and shows the effect. The two plots should be identical except for the scaling of the x-axis, instead they look very different.

import mass
import numpy as np
import pylab as plt
plt.ion()

def plotline(escale, taubigscale):

    fwhm = 2*escale
    tail_tau = 30*escale
    dph_de=1
    def cleanspectrum_fn(x): return mass.getmodel(mass.STANDARD_FEATURES["MnKAlpha"]*escale).spect.pdf(x, instrument_gaussian_fwhm=fwhm)
    bin_width=0.5*escale
    energy = np.arange(5800,6000,0.5)*escale
    length_lo = tail_tau*dph_de/bin_width
    length_hi = 0
    peak_ph = mass.STANDARD_FEATURES["MnKAlpha"]*escale
    tail_frac = 0.5
    tail_frac_hi = 0
    print(f"""{length_lo=} {length_hi=} {fwhm=} {peak_ph=} {dph_de=} {tail_frac=} 
    {tail_tau=} {tail_frac_hi=}  {bin_width=}
    {len(energy)=}""")
    spect = mass.line_fits._smear_exponential_tail(
    cleanspectrum_fn, energy, fwhm, tail_frac, length_lo, tail_frac_hi, length_hi)
    spect_scaletau = mass.line_fits._smear_exponential_tail(
    cleanspectrum_fn, energy, fwhm, tail_frac, length_lo*taubigscale, tail_frac_hi, length_hi)
    plt.figure()
    plt.plot(energy, spect, label=f"tail_tau={tail_tau:.2e}")
    plt.plot(energy, spect_scaletau, label=f"tail_tau={(tail_tau*taubigscale):.2e}")
    plt.legend()
    plt.title(f"{escale=} {taubigscale=}")
    plt.xlabel("energy (eV)")
    plt.ylabel("counts  (arb)")

plotline(1,100)
plotline(1000, 100)

‌

[joefowler]

Original comment by Joseph Fowler (Bitbucket: joe_fowler, ).

---

I can’t run this. Is my mass installation messed up? Are you running the latest?

‌

import mass
import numpy as np
import pylab as plt
plt.ion()

def plotline(escale, taubigscale):

    fwhm = 2*escale
    tail_tau = 30*escale
    dph_de=1
    def cleanspectrum_fn(x): return mass.getmodel(mass.STANDARD_FEATURES["MnKAlpha"]*escale).spect.pdf(x, instrument_gaussian_fwhm=fwhm)
    bin_width=0.5*escale
    energy = np.arange(5800,6000,0.5)*escale
    length_lo = tail_tau*dph_de/bin_width
    length_hi = 0
    peak_ph = mass.STANDARD_FEATURES["MnKAlpha"]*escale
    tail_frac = 0.5
    tail_frac_hi = 0
    print(f"""{length_lo=} {length_hi=} {fwhm=} {peak_ph=} {dph_de=} {tail_frac=} 
    {tail_tau=} {tail_frac_hi=}  {bin_width=}
    {len(energy)=}""")
    spect = mass.line_fits._smear_exponential_tail(
    cleanspectrum_fn, energy, fwhm, tail_frac, length_lo, tail_frac_hi, length_hi)
    spect_scaletau = mass.line_fits._smear_exponential_tail(
    cleanspectrum_fn, energy, fwhm, tail_frac, length_lo*taubigscale, tail_frac_hi, length_hi)
    plt.figure()
    plt.plot(energy, spect, label=f"tail_tau={tail_tau:.2e}")
    plt.plot(energy, spect_scaletau, label=f"tail_tau={(tail_tau*taubigscale):.2e}")
    plt.legend()
    plt.title(f"{escale=} {taubigscale=}")
    plt.xlabel("energy (eV)")
    plt.ylabel("counts  (arb)")

plotline(1,100)
plotline(1000, 100)

## -- End pasted text --
length_lo=60.0 length_hi=0 fwhm=2 peak_ph=5898.801 dph_de=1 tail_frac=0.5 
    tail_tau=30 tail_frac_hi=0  bin_width=0.5
    len(energy)=400
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[1], line 35
     32     plt.xlabel("energy (eV)")
     33     plt.ylabel("counts  (arb)")
---> 35 plotline(1,100)
     36 plotline(1000, 100)

Cell In[1], line 23, in plotline(escale, taubigscale)
     19 tail_frac_hi = 0
     20 print(f"""{length_lo=} {length_hi=} {fwhm=} {peak_ph=} {dph_de=} {tail_frac=} 
     21 {tail_tau=} {tail_frac_hi=}  {bin_width=}
     22 {len(energy)=}""")
---> 23 spect = mass.line_fits._smear_exponential_tail(
     24 cleanspectrum_fn, energy, fwhm, tail_frac, length_lo, tail_frac_hi, length_hi)
     25 spect_scaletau = mass.line_fits._smear_exponential_tail(
     26 cleanspectrum_fn, energy, fwhm, tail_frac, length_lo*taubigscale, tail_frac_hi, length_hi)
     27 plt.figure()

AttributeError: module 'mass' has no attribute 'line_fits'

‌

[joefowler]

Original comment by Joseph Fowler (Bitbucket: joe_fowler, ).

---

Okay, replace line_fits → line_models to fix that.

On the actual bug, I am mystified. Even if the scaling is wrong, why is one orange curve long-tailed and the other one zero-tailed?

I would point out that the exponential smearing is kind of crap. It is done in Fourier space, hence very long tails (tails comparable in size to the energy width of the spectrum) will wrap around. To see this in action, take your two figures made in the example above, and use the “l” key to look at them in semilog-y form. We need the smearing to be fast, so we do it this way.

That's not a refusal-to-fix. Just a side point. I’ll see if I can make sense of this. Clearly some scaling isn’t being used correctly.

[joefowler]

Original comment by Galen O'Neil (Bitbucket: oneilg, GitHub: oneilg).

---

I did all this testing on the TrueBq computer, which has an older mass install. Sorry for the out of dateness.

I did notice some impact on the right side of the plot, which surprised me. Your explanation makes sense.

[joefowler]

Original comment by Joseph Fowler (Bitbucket: joe_fowler, ).

---

Okay, I think the error is in your example, not in Mass. I think you’re interpreting length_lo as a number of bins, but it is supposed to be in energy units (or technically, in the same units as the second (energy) argument to _smear_exponential_tail() So

import mass
import numpy as np
import pylab as plt
plt.ion()

def plotline(escale, taubigscale, tail_frac=0.5):

    fwhm = 2*escale
    tail_tau = 30*escale
    dph_de=1
    def cleanspectrum_fn(x): 
        return mass.getmodel(mass.STANDARD_FEATURES["MnKAlpha"]*escale).spect.pdf(x, instrument_gaussian_fwhm=fwhm)

    bin_width=0.5*escale
    energy = np.arange(5800, 6000, 0.5)*escale
    length_lo = tail_tau*dph_de
    length_hi = 0
    # peak_ph = mass.STANDARD_FEATURES["MnKAlpha"]*escale
    tail_frac_hi = 0

    print(f"""{length_lo=} {length_hi=} {fwhm=} {dph_de=} {tail_frac=} 
        {tail_tau=} {tail_frac_hi=}  {bin_width=}
        {len(energy)=}""")
    spect = mass.line_models._smear_exponential_tail(
        cleanspectrum_fn, energy, fwhm, tail_frac, length_lo, tail_frac_hi, length_hi)
    spect_scaletau = mass.line_models._smear_exponential_tail(
        cleanspectrum_fn, energy, fwhm, tail_frac, length_lo*taubigscale, tail_frac_hi, length_hi)

    plt.figure()
    plt.plot(energy, spect, label=f"tail_tau={tail_tau:.2e}")
    plt.plot(energy, spect_scaletau, label=f"tail_tau={(tail_tau*taubigscale):.2e}")
    plt.legend()
    plt.title(f"{escale=} {taubigscale=}")
    plt.xlabel("energy (eV)")
    plt.ylabel("counts  (arb)")

plotline(1,100)
plotline(1000, 100)

If you change your example function to the above version, I think your two figures come out the same (except they differ by x1000 in the x-axis).

So I don’t agree that your latest MRE is a bug. However, the previous one (“Ok I made a better minimum reproducing example…”) is still troubling.

[joefowler]

Original comment by Galen O'Neil (Bitbucket: oneilg, GitHub: oneilg).

---

In that case I think the bug is here: https://bitbucket.org/joe_fowler/mass/src/253a07a1fc290191b194948f99a2b188437662f3/mass/calibration/line_models.py#lines-229

‌

237 and 238, which is where I drew the definition of length_lo and length_hi from. I’m pretty sure I wrote that code, but copied most of it from the previous linefit code. Do I just remove dividing out bin width?

[joefowler]

Original comment by Galen O'Neil (Bitbucket: oneilg, GitHub: oneilg).

---

fix issue 250 and update tests

[joefowler]

Original comment by Galen O'Neil (Bitbucket: oneilg, GitHub: oneilg).

---

wtf bitbucket, me opening a pull request does not resolve this issue, I haven’t merged the pull request!!