-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathelastic_ace_endl_native.py
More file actions
195 lines (152 loc) · 8.11 KB
/
elastic_ace_endl_native.py
File metadata and controls
195 lines (152 loc) · 8.11 KB
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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
#! /usr/bin/env python
import PyFrensie.Data.Native as Native
import PyFrensie.Data.ACE as ACE
import PyFrensie.Data.ENDL as ENDL
import PyFrensie.Utility as Utility
import PyFrensie.Utility.Prng as Prng
import PyFrensie.Utility.Interpolation as Interpolation
import PyFrensie.MonteCarlo.Collision as Collision
import PyTrilinos.Teuchos as Teuchos
import numpy
import matplotlib.pyplot as plt
datadir = '/home/software/mcnpdata/'
#datadir = '/home/lkersting/frensie/src/packages/test_files/'
source = Teuchos.FileInputSource( datadir + '/cross_sections.xml' )
xml_obj = source.getObject()
cs_list = Teuchos.XMLParameterListReader().toParameterList( xml_obj )
energy = 15.7
# -------------------------------------------------------------------------- ##
# Native Data
# -------------------------------------------------------------------------- ##
data_list = cs_list.get( 'Au-Native' )
file_name = datadir + data_list.get( 'electroatomic_file_path' )
native_data = Native.ElectronPhotonRelaxationDataContainer( file_name )
native_energy_grid = native_data.getElectronEnergyGrid()
analog_dist = Collision.createAnalogElasticDistribution(native_data, True, True, 1e-15)
#lin_analog_dist = Collision.createAnalogElasticDistribution(native_data, False, True, 1e-15)
native_index = 0
for i in range(0, native_energy_grid.size ):
if native_energy_grid[i] <= energy:
native_index = i
energy_0 = native_energy_grid[native_index]
energy_1 = native_energy_grid[native_index+1]
native_lin_interp = ( energy - energy_0 )/( energy_1 - energy_0 )
native_log_interp = numpy.log(energy/energy_0)/numpy.log(energy_1/energy_0)
# -------------------------------------------------------------------------- ##
# Native Cutoff Elastic Cross Section
# -------------------------------------------------------------------------- ##
elastic_cross_sections = native_data.getCutoffElasticCrossSection()
cs_0 = elastic_cross_sections[native_index]
cs_1 = elastic_cross_sections[native_index+1]
native_elastic_cs = cs_0 + (cs_1 - cs_0)*native_lin_interp
#native_elastic_cs = numpy.exp(numpy.log(cs_0) + numpy.log(cs_1/cs_0)*native_log_interp)
# -------------------------------------------------------------------------- ##
# Native Analog Elastic Cross Section
# -------------------------------------------------------------------------- ##
native_analog_cs = native_elastic_cs/analog_dist.evaluateCDFAtCutoff( energy )
#native_elastic_cs = native_elastic_cs/lin_analog_dist.evaluateCDFAtCutoff( energy )
# -------------------------------------------------------------------------- ##
# ACE Data
# -------------------------------------------------------------------------- ##
data_list = cs_list.get( 'Au' )
file_name = datadir + data_list.get( 'electroatomic_file_path' )
table_name = data_list.get( 'electroatomic_table_name' )
table_start = data_list.get( 'electroatomic_file_start_line' )
ace_file = ACE.ACEFileHandler( file_name, table_name, table_start )
xss_extractor = ACE.XSSEPRDataExtractor( ace_file.getTableNXSArray(), ace_file.getTableJXSArray(), ace_file.getTableXSSArray() )
ace_energy_grid = xss_extractor.extractElectronEnergyGrid()
ace_index = 0
for i in range(0, ace_energy_grid.size ):
if ace_energy_grid[i] <= energy:
ace_index = i
energy_0 = ace_energy_grid[ace_index]
energy_1 = ace_energy_grid[ace_index+1]
lin_interp = ( energy - energy_0 )/( energy_1 - energy_0 )
log_interp = numpy.log(energy/energy_0)/numpy.log(energy_1/energy_0)
# -------------------------------------------------------------------------- ##
# ACE Cutoff Elastic Cross Section
# -------------------------------------------------------------------------- ##
elastic_cross_sections = xss_extractor.extractElasticCrossSection()
cs_0 = elastic_cross_sections[ace_index]
cs_1 = elastic_cross_sections[ace_index+1]
#ace_elastic_cs = cs_0 + (cs_1 - cs_0)*lin_interp
ace_elastic_cs = numpy.exp(numpy.log(cs_0) + numpy.log(cs_1/cs_0)*log_interp)
# -------------------------------------------------------------------------- ##
# ENDL Data
# -------------------------------------------------------------------------- ##
data_list = cs_list.get( 'Au-ENDL' )
file_name = datadir + data_list.get( 'electroatomic_file_path' )
endl_data = ENDL.ENDLDataContainer( file_name )
# -------------------------------------------------------------------------- ##
# ENDL Cutoff Elastic Cross Section
# -------------------------------------------------------------------------- ##
elastic_cs = endl_data.getCutoffElasticCrossSection()
endl_energy_grid = endl_data.getElasticEnergyGrid()
endl_index = 0
for n in range(0, endl_energy_grid.size ):
if endl_energy_grid[n] <= energy:
endl_index = n
energy_0 = endl_energy_grid[endl_index]
energy_1 = endl_energy_grid[endl_index+1]
lin_interp = ( energy - energy_0 )/( energy_1 - energy_0 )
log_interp = numpy.log(energy/energy_0)/numpy.log(energy_1/energy_0)
cs_0 = elastic_cs[endl_index]
cs_1 = elastic_cs[endl_index+1]
#endl_elastic_cs = cs_0 + (cs_1 - cs_0)*lin_interp
endl_elastic_cs = numpy.exp(numpy.log(cs_0) + numpy.log(cs_1/cs_0)*log_interp)
# -------------------------------------------------------------------------- ##
# ENDL Analog Cross Section
# -------------------------------------------------------------------------- ##
analog_cs = endl_data.getTotalElasticCrossSection()
cs_0 = analog_cs[endl_index]
cs_1 = analog_cs[endl_index+1]
#endl_analog_cs = cs_0 + (cs_1 - cs_0)*lin_interp
endl_analog_cs = numpy.exp(numpy.log(cs_0) + numpy.log(cs_1/cs_0)*log_interp)
angular_energy_grid = endl_data.getCutoffElasticAngularEnergyGrid()
angular_index = 0
for n in range(0, angular_energy_grid.size ):
if angular_energy_grid[n] <= energy:
angular_index = n
energy_0 = angular_energy_grid[angular_index]
energy_1 = angular_energy_grid[angular_index+1]
lin_interp = ( energy - energy_0 )/( energy_1 - energy_0 )
log_interp = numpy.log(energy/energy_0)/numpy.log(energy_1/energy_0)
pdf_0 = endl_data.getCutoffElasticPDFAtEnergy(energy_0)
cutoff_pdf_0 = pdf_0[0]
pdf_1 = endl_data.getCutoffElasticPDFAtEnergy(energy_1)
cutoff_pdf_1 = pdf_1[0]
cutoff_pdf = cutoff_pdf_0 + (cutoff_pdf_1 - cutoff_pdf_0)*log_interp
cutoff_pdf_lin = cutoff_pdf_0 + (cutoff_pdf_1 - cutoff_pdf_0)*lin_interp
cutoff_pdf_loglog = numpy.exp(numpy.log(cutoff_pdf_0) + numpy.log(cutoff_pdf_1/cutoff_pdf_0)*log_interp)
eta = analog_dist.evaluateMoliereScreeningConstant( energy )
cutoff_cdf = eta/( eta + cutoff_pdf*(eta*1.0e-6 + 1.0e-12) )
cutoff_cdf_lin = eta/( eta + cutoff_pdf*(eta*1.0e-6 + 1.0e-12) )
cutoff_cdf_loglog = eta/( eta + cutoff_pdf*(eta*1.0e-6 + 1.0e-12) )
endl_analog_cs_2 = endl_elastic_cs/cutoff_cdf
endl_analog_cs_3 = endl_elastic_cs/cutoff_cdf_lin
endl_analog_cs_4 = endl_elastic_cs/cutoff_cdf_loglog
# -------------------------------------------------------------------------- ##
# Elastic Data Comparison
# -------------------------------------------------------------------------- ##
print "\n----- Cutoff Elastic Cross Section at",energy,"-----"
print '\tNative cs = ','%.16e' % native_elastic_cs
print '\tACE cs = ','%.16e' % ace_elastic_cs
print '\tENDL cs = ','%.16e' % endl_elastic_cs
rel_diff = (native_elastic_cs-endl_elastic_cs)/endl_elastic_cs
print '\tRel diff (Native vs ENDL) = ','%.16e' % rel_diff
rel_diff = (ace_elastic_cs-endl_elastic_cs)/endl_elastic_cs
print '\tRel diff (ACE vs ENDL) = ','%.16e' % rel_diff
print "\n----- Analog Elastic Cross Section at",energy,"-----"
print '\tNative cs = ','%.16e' % native_analog_cs
print '\tENDL cs (tables) = ','%.16e' % endl_analog_cs
print '\tENDL cs (calc) = ','%.16e' % endl_analog_cs_2
print '\tENDL cs (calc) = ','%.16e' % endl_analog_cs_3
print '\tENDL cs (calc) = ','%.16e' % endl_analog_cs_4
rel_diff = (native_analog_cs-endl_analog_cs)/endl_analog_cs
print '\trel diff (Native vs ENDL Tables) = ','%.16e' % rel_diff
rel_diff = (endl_analog_cs_2-endl_analog_cs)/endl_analog_cs
print '\trel diff (lin-log ENDL Calc vs ENDL Tables) = ','%.16e' % rel_diff
rel_diff = (endl_analog_cs_3-endl_analog_cs)/endl_analog_cs
print '\trel diff (lin-lin ENDL Calc vs ENDL Tables) = ','%.16e' % rel_diff
rel_diff = (endl_analog_cs_4-endl_analog_cs)/endl_analog_cs
print '\trel diff (log-log ENDL Calc vs ENDL Tables) = ','%.16e' % rel_diff