#!/usr/bin/env python3
#
'''
.. module:: MCeventlist
:synopsis: make Monte-Carlo event list for Holzer alpha1/2 and/or Beta lines spectrum
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import numpy
from tesfdmtools.methods.holzer import holzer,hbeta,holzerelements
from tesfdmtools.methods.brems import brems
#from holzer import holzer,hbeta,holzerelements
import matplotlib.pyplot as plt
def plotspec(x,y):
'''
show spectrum
Args:
* `x` = energy axis (eV)
* `y` = counts
'''
f,ax=plt.subplots(1)
ax.plot(x,y)
ax.set_xlabel('Energy (eV)')
ax.set_ylabel('Counts')
ax.set_title('Generated spectrum')
plt.show()
plt.close('all')
[docs]def MCeventlist(erange=None,counts=2000,res=2.5,elements=['mn'],lines=['Alpha'],bspec=None):
'''
make Monte-Carlo event list for Holzer alpha1/2 or beta lines spectrum
Kwargs:
* `erange` = range in eV for which the eventlist should be made
* `counts` = number of events to generate
* `res` = resolution of spectrum in eV
* `elements` = list of elements [ti,cr,mn,cu] for which to show lines.
* `lines` = list of line types [Alpha,Beta] to show
* `bspec` = add Bremsstralung spectrum model with given factor if not **None**
Returns:
* `elist` = list of event energies
* `spectrum` = input spectrum for the events
* `eax` = energy axis for the spectrum
'''
if ( lines is None ) or ( len(lines) == 0 ):
lines=['Alpha']
if ( erange is None ) or ( erange[0] != erange[1] ):
erangeh=holzerelements(elements=elements)
emin=9E9
emax=-9E9
for line in erangeh:
emin=min([emin,min(erangeh[line])])
emax=max([emax,max(erangeh[line])])
if erange is None:
erange=[emin-40.0,emax+40.0]
else:
erange=[5850.0,5940.0]
bins=int((erange[1]-erange[0])/(res/500.0))
bins=2**int(numpy.log(bins)/numpy.log(2.0)+1.0) # make it power of two for fast fft later
binnr=numpy.arange(bins,dtype=float)
eax=binnr/float(bins)*(erange[1]-erange[0])+erange[0]
hpar=numpy.array([0.0,1.0,1000.0,(res/1000.0)],dtype=float)
for j,line in enumerate(lines):
if line == 'Alpha':
if j == 0:
spectrum=holzer(hpar,(eax/1000.0))
else:
spectrum=spectrum+holzer(hpar,(eax/1000.0))
elif line == 'Beta':
if j == 0:
spectrum=0.2*hbeta(hpar,(eax/1000.0))
else:
spectrum=spectrum+0.2*hbeta(hpar,(eax/1000.0))
if bspec is not None:
if isinstance(bspec,bool):
bspec=0.709 # default: ratio between MXS line intensities and Brems. spectrum
spectrum=spectrum+brems(eax)*bspec*numpy.sum(spectrum)
# plotspec(eax,spectrum) # for debugging purposes, show spectrum
cumsp=numpy.cumsum(spectrum) # cumulative spectrum
cumsp=cumsp/cumsp[-1]*float(bins)
invsp=numpy.interp(binnr,cumsp,eax) # invert cumulative spectrum
elist=invsp[numpy.array((numpy.random.random(counts)*bins),dtype=int)] # Monte-Carlo
return (elist,spectrum,eax) # return list of events and input spectrum
#==========================================================================================
if __name__ == "__main__":
'''
Test program
'''
import sys
# import matplotlib.pyplot as plt
elements=['mn']
lines=['Alpha']
bspec=None
if len(sys.argv) > 1:
elements=[]
lines=[]
for pp in sys.argv[1:]:
if pp in ['Alpha','Beta']:
lines.append(pp)
elif pp in ['mn','cr','cu','ti']:
elements.append(pp)
elif pp == 'bspec':
bspec=True
else:
sys.exit("Use: %s element1 [ element2 ....] line1 [ line2 ]" % sys.argv[0] )
if len(lines) == 0:
lines=['Alpha']
if len(elements) == 0:
elements=['mn']
elist,ss,ee=MCeventlist(counts=10000,elements=elements,lines=lines,bspec=bspec)
f,ax=plt.subplots(1)
ax.set_xlabel('Energy (eV)')
ax.set_ylabel('Counts')
ax.set_title("Spectrum: '"+"' '".join(elements)+"'")
ax.hist(elist,bins=2048)
plt.show()
plt.close('all')