/usr/lib/python2.7/dist-packages/PyMca/EPDL97/GenerateEPDL97TotalCrossSections.py is in pymca 4.7.1+dfsg-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 | __doc__= "Generate specfile from EPL97 total cross sections in keV and barn"
import os
import sys
import EPDL97Parser as EPDLParser
Elements = EPDLParser.Elements
AVOGADRO_NUMBER = EPDLParser.AVOGADRO_NUMBER
import numpy
log = numpy.log
exp = numpy.exp
getTotalCoherentCrossSection = EPDLParser.getTotalCoherentCrossSection
getTotalIncoherentCrossSection = EPDLParser.getTotalIncoherentCrossSection
getTotalPhotoelectricCrossSection = EPDLParser.getTotalPhotoelectricCrossSection
getTotalPairCrossSection = EPDLParser.getTotalPairCrossSection
getTotalTripletCrossSection = EPDLParser.getTotalTripletCrossSection
if len(sys.argv) < 3:
print("Usage:")
print("python EPDLGenerateTotalCrossSections SPEC_output_filename barns_flag")
sys.exit(0)
def getHeader(filename):
text = '#F %s\n' % filename
text += '#U00 This file is a direct conversion to specfile format of \n'
text += '#U01 the original EPDL97 total cross sections contained in the\n'
text += '#U02 EPDL97.DAT from the library.\n'
text += '#U03 EPDL97 itself can be found at:\n'
text += '#U04 http://www-nds.iaea.org/epdl97/libsall.htm\n'
text += '\n'
return text
fname = sys.argv[1]
if os.path.exists(fname):
os.remove(fname)
if int(sys.argv[2]):
BARNS = True
else:
BARNS = False
print("BARNS = %s" % BARNS)
outfile = open(fname, 'wb')
outfile.write(getHeader(fname))
for i in range(1, 101):
print("i = %d element = %s" % (i, Elements[i-1]))
#coherent
energy_cohe, value_cohe, mode_cohe = getTotalCoherentCrossSection(i,
getmode=True)
#incoherent
energy_incohe, value_incohe, mode_incohe = getTotalIncoherentCrossSection(i,
getmode=True)
#photoelectric
energy_photo, value_photo, mode_photo = getTotalPhotoelectricCrossSection(i,
getmode=True)
#check to see the energies:
#for j in range(10):
# print energy_cohe[j], energy_incohe[j], energy_photo[j]
#to select an appropriate energy grid as close as possible to the original
#while keeping in mind the PyMca goals, I use the coherent energy grid till
#the non-zero first value of the photoelectric cross section. At that point,
#I use the photoelectric energy grid.
energy = numpy.concatenate((energy_cohe[energy_cohe<energy_photo[0]],
energy_photo))
#now perform a log-log interpolation when needed
#lin-lin interpolation:
#
# y0 (x1-x) + y1 (x-x0)
# y = -------------------------
# x1 - x0
#
#log-log interpolation:
#
# log(y0) * log(x1/x) + log(y1) * log(x/x0)
# log(y) = ------------------------------------------
# log (x1/x0)
#
cohe = numpy.zeros(len(energy), numpy.float)
incohe = numpy.zeros(len(energy), numpy.float)
photo = numpy.zeros(len(energy), numpy.float)
total = numpy.zeros(len(energy), numpy.float)
#coherent needs to interpolate
indices = numpy.nonzero(energy_cohe<energy_photo[0])
cohe[indices] = value_cohe[indices]
for n in range(len(indices),len(energy)):
x = energy[n]
j1 = len(indices)
while energy_cohe[j1] < x:
j1 += 1
j0 = j1 - 1
x0 = energy_cohe[j0]
x1 = energy_cohe[j1]
y0 = value_cohe[j0]
y1 = value_cohe[j1]
cohe[n] = exp((log(y0) * log(x1/x) + log(y1) * log(x/x0))/log(x1/x0))
#compton needs to interpolate everything
for n in range(len(energy)):
x = energy[n]
j1 = 0
while energy_incohe[j1] < x:
j1 += 1
j0 = j1 - 1
x0 = energy_incohe[j0]
x1 = energy_incohe[j1]
y0 = value_incohe[j0]
y1 = value_incohe[j1]
incohe[n] = exp((log(y0) * log(x1/x) + log(y1) * log(x/x0))/log(x1/x0))
#photoelectric does not need to interpolate anything
photo[energy>=energy_photo[0]] = value_photo[:]
#convert to keV and cut at 500 keV
energy *= 1000.
indices = numpy.nonzero(energy<=500.)
energy = energy[indices]
photo = photo[indices]
cohe = cohe[indices]
incohe = incohe[indices]
#I cut at 500 keV, I do not need to take the pair production
total = photo + cohe + incohe
#now I am ready to write a Specfile
ele = Elements[i-1]
text = '#S %d %s\n' % (i, ele)
text += '#N 5\n'
labels = '#L PhotonEnergy[keV]'
labels += ' Rayleigh(coherent)[barn/atom]'
labels += ' Compton(incoherent)[barn/atom]'
labels += ' CoherentPlusIncoherent[barn/atom]'
labels += ' Photoelectric[barn/atom]'
labels += ' TotalCrossSection[barn/atom]\n'
if not BARNS:
labels = labels.replace("barn/atom", "cm2/g")
factor = (1.0E-24*AVOGADRO_NUMBER)/EPDLParser.getAtomicWeights()[i-1]
else:
factor = 1.0
text += labels
if 0:
fformat = "%g %g %g %g %g %g\n"
else:
fformat = "%.6E %.6E %.6E %.6E %.6E %.6E\n"
outfile.write(text)
for n in range(len(energy)):
line = fformat % (energy[n],
cohe[n] * factor,
incohe[n] * factor,
(cohe[n]+incohe[n]) * factor,
photo[n] * factor,
total[n] * factor)
outfile.write(line)
outfile.write('\n')
outfile.close()
|