Coverage for larch/xsw/fluo_det.py: 0%
450 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
1"""
2Fluorescence Intensity calculations:
3 Y. Choi
5calculates fluorescence intensities for multiple elements at a fixed incident x-ray energy.
6measured fluorescence depends on net fluorescence yield and transmission.
7included are:
8 (with incident energy dependence): atomic crosssection, fluorescence yield for each edge,
9 transition probability for each emission line,
10 (with emitted fluorescence energy dependence): transmission to vortex detector, detector efficiency.
11global variables: re, Barn, Nav, pre_edge_margin, fluo_emit_min, det_res
12AtNum2f1f2Xect: needs datafile readf1f2a.py which is from chantler table
13cal_NetYield2: for XRF calculations. major and minor lines are calculated
14 output as 'Elemental_Sensitivity.txt'
16sim_spectra: keeps individual emission lines with significant intensity without weight-averaging
17 each element can have a different relative concentration.
18 makes strongest emission to 100 and scales the rest accordingly.
19 output as 'simSpectrum.txt'
20 runs sim_GaussPeaks() to emulate measure intensity as gaussian peak
21 output as 'simSpectrum.plot'
22 runs cal_NetYield2()
235/24/2010: self-absorption effect included, 4 options (top, bottom, all, surface)
24 Each option makes different elemental distribution,
25 but substrate elements are always distributed evenly.
26 get_ChemName(), nominal_density() are in readf1f2a.py
275/26/2010: incident x-ray attenuation effect added, and incident angle effect is also added.
28 refraction effect added to the angle, below critical angle, th=1e-10.
296/4/2010: SampleMatrix2 has two layers and each layer can have different composition
30 Class ElemFY is for fluorescence elements with AtSym, Conc, tag attribues,
31 inputs in SimSpect, Cal_NetFyield2 changed.
32 Concentrations are normalized to the concentration of the substrate1 molecule.
33 Ex) with CaCO3[2.71g/cc]/Fe2O3[5.26g/cc] substrate, Ca concentration is 1.0 meaning 1.63e22 Ca/cc
34 and Fe concentration is 2.43 = 2*1.98e22/1.63e22 where 1.98e22 is number density (/cc) of Fe2O3
356/8/2010: now using python2.6 instead of 2.5
364/14/2011: fluo_det.py: detection dependent part in fluorescence yield.
37 Original fluo.py was modified to fluo_new.py.
38 readf1f2a.py is used instead of original readf1f2.py
39 fluo_new.py is being split into two.
40 fluo_elem.py is for elemtnal dependence: fluorescenc yield, crosssection.
41 fluo_det.py is for detection dependnece: detector efficiency, attenuation, sample
42 * for attenuation, total cross-section([4]) is used instead of photoelectric ([2]).
43 ** for fluorescence yield, photoelectric is used.
44"""
47import math
48import numpy
49import sys
50from larch.utils.physical_constants import AVOGADRO, BARN
51from larch.xray import xray_delta_beta, chemparse, atomic_mass
53pre_edge_margin=150. # FY calculated from 150 eV below the absorption edge.
54fluo_emit_min=500. # minimum energy for emitted fluorescence. ignore fluorescence emissions below 500eV
55det_res=100. # detector resoltuion in eV, used in sim_GaussPeaks
57'''
58----------------------------------------------------------------------------------------------------------
59class: Material, ElemFY, SampleMatrix2
60----------------------------------------------------------------------------------------------------------
61'''
62class Material:
63 def __init__(self, composition, density, thickness=1):
64 self.composition=composition # ex) SiO2
65 self.density=density # in g/cm^3
66 self.thickness=thickness # in cm
67 self.la=1.0 #absoprtion length in cm
68 self.trans=0.5 # transmission.
69 self.absrp=0.5 # absorption
70 self.delta=0.1 # index of refraction, real part
71 self.beta=0.1 # index of refraction, imaginary part
72 # MN replace
73 elements = chemparse(composition)
74 AtomWeight=0
75 # AtomList = elements.keys()
76 # AtomIndex = elements.values()
77 # for (ii, atom) in enumerate(AtomList):
78 # # MN replace...
79 # AtWt= atomic_mass(atom)
80 # index=AtomIndex[ii]
81 # AtomWeight = AtomWeight + index*AtWt
82 for sym, mole in elements.items():
83 AtomWeight = AtomWeight + mole*atomic_mass(sym)
84 NumberDensity = density*AVOGADRO/AtomWeight
85 self.NumDen = NumberDensity # number of molecules per cm^3
86 self.AtWt = AtomWeight # weight per mole
88 def getLa(self, energy, NumLayer=1.0): # get absorption legnth for incident x-ray, NumLayer for multiple layers
89 # MN replace...
90 temp= xray_delta_beta(self.composition, self.density, energy)
91 # returns delta, beta, la_photoE, Nat, la_total
92 self.delta=temp[0]; self.beta=temp[1]
93 self.la=temp[2] # temp[4] (total) instead of temp[2] (photoelectric)
94 if NumLayer<0: NumLayer=0 # Number of layers cannot be less than zero
95 self.trans=math.exp(-self.thickness*NumLayer/self.la) #la attenuation length cm, NumLayer for multiple layers
96 self.absrp=1-self.trans
99class ElemFY: # fluorescing element
100 def __init__(self, AtomicSymbol='Fe', Concentration=10e-6, tag='dilute'):
101 temp=AtomicSymbol[0].upper() # capitalize the first letter
102 self.AtSym=temp+AtomicSymbol[1:]
103 self.Conc=Concentration
104 self.tag=tag # dilute, substrate1, substrate2
107class SampleMatrix2: # sample matrix for self-absorption correction, 6/3: two layers with different compositions
108 def __init__(self, composition1='Si', density1=2.33, thickness1=0.1, \
109 composition2='Si', density2=2.33, thickness2=0.1,
110 angle0=45.0, option='surface'):
111 self.composition1 = composition1 # ex) Fe2O3
112 # MN replace:
113 out= chemparse(composition1) # output from get_ChemName
114 self.ElemList1 = out.keys() # list of elments in matrix: 'Fe', 'O' for Fe2O3
115 self.ElemInd1 = out.values() # list of index: 2, 3 for Fe2O3
116 # self.ElemFrt1 = out[2] # list of fraction: 0.4, 0.6 for Fe2O3
117 self.density1 = density1 # in g/cm^3
118 self.thickness1 = thickness1 # top layer thickness in cm
119 self.composition2 = composition2
120 # MN replace with chemparse()
121 out=chemparse(composition2)
122 self.ElemList2 = out.keys() # for the bottom substrate
123 self.ElemInd2 = out.values()
124 # self.ElemFrt2 = out[2]
125 self.density2 = density2
126 self.thickness2 = thickness2 # bottom layer thickness in cm
127 self.angle = angle0*math.pi/180. # in radian, incident beam angle, surface normal =pi/2
128 self.option = option # option for fluorescing element location, surface/top/bottom
129 self.la1 = 1.0 # absoprtion length in cm
130 self.delta1 = 0.1 # index of refraction, real part correction
131 self.beta1 = 0.1 # index of refraction, imaginary part
132 self.Nat1 = 1.0 # atomic number density of the Fe2O3, 1e22 Fe2O3 atoms /cm^3
133 self.la2 = 1.0
134 self.delta2 = 0.1
135 self.beta2 = 0.1
136 self.Nat2 = 1.0 # atomic number density of the first element, 1e22 Fe2O3 atoms/cm^3
137 self.scale = 1.0 # weighted average over the depth range: sum of (trans*factors)
138 self.scale1 = 1.0 # sum of (trans*1.0) this is for substrate element in top layer
139 self.scale2 = 1.0 # sum of (trans*thickness2/thickness1) for substrate element in bottom layer
140 self.ElemListFY =[] # fluorescing element
141 self.ElemFrtFY =[] # fraction for fluorescing element
142 self.Nat1=1 # atomic number density in atoms/cc, for example # of Fe2O3/cm^3 for Fe2O3 substrate
143 self.Nat2=1 # atomic number density in atoms/cc
144 substrate1_material = Material(composition1, density1)
145 AtNumDen1 = substrate1_material.NumDen # substrate1
146 substrate2_material = Material(composition2, density2)
147 AtNumDen2 = substrate2_material.NumDen # substrate2, fixed 8/12/10
148 self.Nat1=AtNumDen1
149 self.Nat2=AtNumDen2
150 self.txt=''
151 text1='substrate1:%6.3e %s/cm^3' % (AtNumDen1, composition1)
152 #print(text1)
153 text2='substrate2:%6.3e %s/cm^3' % (AtNumDen2, composition2)
154 self.txt=text1+' '+text2
155 print(self.txt)
156 # atom.conc is normalized to the number density of substrate1 molecule
157 for (ii, item) in enumerate(self.ElemList1):
158 # MN replace:
159 if atomic_number(item)>=12: # ignore elements below Mg
160 #atom=ElemFY(item, self.ElemFrt1[ii], 'substrate1')
161 atom=ElemFY(item, self.ElemInd1[ii], 'substrate1')
162 # eg. for Fe2O3, Fe concentration is 2 (=2 x Fe2O3 number density)
163 self.ElemListFY.append(atom)
164 for (ii, item) in enumerate(self.ElemList2):
165 # MN replace:
166 if atomic_number(item)>=12: # ignore elements below Mg
167 #atom=ElemFY(item, self.ElemFrt2[ii], 'substrate2')
168 atom=ElemFY(item, self.ElemInd2[ii]*AtNumDen2/AtNumDen1, 'substrate2')
169 self.ElemListFY.append(atom)
170 numLayer=100 # each layer is sliced into 100 sublayers.
171 self.depths=[] # depth values for calculation
172 for ii in range(numLayer):
173 step=1.0/numLayer
174 self.depths.append(step*ii*self.thickness1)
175 for ii in range(numLayer):
176 step=1.0/numLayer
177 self.depths.append(step*ii*self.thickness2+self.thickness1)
178 self.factors=[] # 1 or pre if element present at each depth
179 pre=self.thickness2/self.thickness1 # prefactor, if two layers have different thickness
180 if self.option=='surface': # fluorescecing atoms present in only on the surface
181 for ii in range(len(self.depths)):
182 if ii==0:
183 self.factors.append(1.0)
184 else:
185 self.factors.append(0.0)
186 if self.option=='all': # fluorescecing atoms present throughout
187 for ii in range(len(self.depths)):
188 if self.depths[ii]<self.thickness1:
189 self.factors.append(1.0)
190 else:
191 self.factors.append(pre)
192 if self.option=='top': # fluorescecing atoms present only in the top layer
193 for ii in range(len(self.depths)):
194 if self.depths[ii]<self.thickness1:
195 self.factors.append(1.0)
196 else:
197 self.factors.append(0.0)
198 if self.option=='bottom': # fluorescecing atoms present only in the bottom layer
199 for ii in range(len(self.depths)):
200 if self.depths[ii]<self.thickness1:
201 self.factors.append(0.0)
202 else:
203 self.factors.append(pre)
204 # note: fluorescing substrate atoms present throughout regardless of the option.
205 self.trans=[] # transmission to surface for emitted fluorescence
206 self.absrp=[] # absorption until surface for emitted fluorescence
207 self.inten0=[] # incident x-ray intensity at each depth
208 for ii in range(len(self.depths)):
209 self.trans.append(0.5)
210 self.absrp.append(0.5)
211 self.inten0.append(0.5)
212 def getPenetration(self, energy0): # incident x-ray penetration(attenuation)
213 # refraction at air/top layer interface
214 # MN replace:
215 temp=f1f2.get_delta(self.composition1, self.density1, energy0) # energy0: incident x-ray energy
216 delta1=temp[0]; beta1=temp[1]
217 la1=temp[4] # in cm, temp[4] instead of temp[2], using total instead of photoelectric
218 self.la1=la1 # absorption length in microns at incident x-ray energy
219 angle_critical1 = (2.0*delta1)**(0.5) # in radian, critical angle for total external reflection
220 if angle_critical1>=self.angle: # below critical angle, the corrected should be zero.
221 angle_corrected1=1.0e-15 # a smaller number instead of zero
222 else: # above critical angle
223 angle_corrected1 = (self.angle**2.0 - angle_critical1**2.0)**(0.5) # in radian
224 # refraction at top/bottom layers interface
225 # MN replace:
226 temp=f1f2.get_delta(self.composition2, self.density2, energy0) # energy0: incident x-ray energy
227 delta2=temp[0]; beta2=temp[1];
228 la2=temp[4] # in cm, temp[4] instead of temp[2], using total instead of photoelectric
229 self.la2=la2 # absorption length in cm at incident x-ray energy
230 angle_corrected2 = ( 2.0-(1.0-delta1)/(1.0-delta2)*(2.0-angle_corrected1**2) )**0.5
231 # using Snell's law, assume beta effect not ignificant, in radian
232 for (ii, depth) in enumerate(self.depths):
233 if self.depths[ii]<self.thickness1: # top layer
234 beampath = depth/math.sin(angle_corrected1) # in cm
235 inten0 = math.exp(-beampath/la1) # attenuated incident beam intensity at depths
236 else: # bottom layer
237 beampath1 = self.thickness1/math.sin(angle_corrected1)
238 beampath2 = (depth-self.thickness1)/math.sin(angle_corrected2)
239 inten0 = math.exp(-beampath1/la1)*math.exp(-beampath2/la2)
240 self.inten0[ii] = inten0 # incident x-ray attenuation
241 def getLa(self, energy, NumLayer=1.0): # emitted fluorescence trasmission attenuation up to top surface
242 transmitted=1.0
243 # MN replace:
244 temp=f1f2.get_delta(self.composition1, self.density1, energy) # energy is for fluorescence
245 self.delta1 = temp[0]
246 self.beta1 = temp[1]
247 self.la1 = temp[4] # in cm, temp[4] instead of temp[2], using total instead of photoelectric
248 # absorption length in cm at emitted fluorescence energy
249 # temp[3]: atomic number density of the first element in atoms/cc
250 # MN replace:
251 temp=f1f2.get_delta(self.composition2, self.density2, energy) # energy is for fluorescence
252 self.delta2 = temp[0]
253 self.beta2 = temp[1]
254 self.la2 = temp[4] # absorption length in cm at emitted fluorescence energy
255 # in cm, temp[4] instead of temp[2], using total instead of photoelectric
256 angle_exit = (math.pi/2.0 - self.angle) # becomes 90 at a small incident angle
257 for (ii, depth) in enumerate(self.depths):
258 if self.depths[ii]<self.thickness1: # top layer
259 transmitted=math.exp(-depth/math.sin(angle_exit)*NumLayer/self.la1)
260 else: # bottom layer
261 transmitted2 = math.exp(-(depth-self.thickness1)/math.sin(angle_exit)*NumLayer/self.la2)
262 transmitted1 = math.exp(-self.thickness1/math.sin(angle_exit)*NumLayer/self.la1)
263 transmitted = transmitted2*transmitted1
264 self.trans[ii] = transmitted
265 self.absrp[ii] = 1.0 - transmitted
266 scale = 0.0; scale1=0.0; scale2=0.0
267 for (ii,trans) in enumerate(self.trans):
268 scale = scale + trans*self.inten0[ii]*self.factors[ii]
269 if self.depths[ii]<self.thickness1:
270 scale1 = scale1 + trans*self.inten0[ii]*1.0
271 else: # if thickness2 is different from thickness1, weight differently
272 scale2 = scale2 + trans*self.inten0[ii]*(self.thickness2/self.thickness1)
273 # scale, scale1, scale2: emitted fluorescence transmission and depth profile
274 self.scale = scale # sum of (trans*factors) for nonsubstrate FY
275 self.scale1 = scale1 # sum of (trans*factors) for substrate FY. factors=1, top layer
276 self.scale2 = scale2 # sum of (trans*factors) for substrate FY. factors=pre, bttom layer
279'''
280----------------------------------------------------------------------------------------------------------
281Detector efficiency, fluorescence attenuation
282----------------------------------------------------------------------------------------------------------
283'''
284# WD30 used for XSW, sample-->He-->Kapton-->Collimator-->Detector
285# WD60 used for XRM, sample-->Collimator-->Detector
286# eV1 is fluorescence energy in eV
287# xHe=1 means Helium gas from sample to WD60 collimator.
288# xAl=1 means 1 layer of 1.5mil Al foil as attenuator
289# xKapton=1 means 1 layer of 0.3 mil Kapton (
290# WD=6 means working distance of 6cm for the detector.
292def Assemble_QuadVortex(eV1):
293 # quad vortex detector efficiency. eV1: fluo energy
294 net=1.
295 BeVortex=Material('Be', 1.85, 0.00125)
296 SiO2Vortex=Material('SiO2', 2.2, 0.00001)
297 SiVortex=Material('Si', 2.33, 0.035)
298 BeVortex.getLa(eV1, 1) # one Be layer in Vortex
299 SiO2Vortex.getLa(eV1, 1) # oxide layer on Si detection layer
300 SiVortex.getLa(eV1, 1) # Si detection layer, what's absorbed is counted.
301 net=net*BeVortex.trans*SiO2Vortex.trans*SiVortex.absrp
302 if (print2screen):
303 print( '%.3f eV : BeVortex.trans=%.3e , SiO2Vortex.trans=%.3e, SiVortex.absrp=%.3e, Det_efficiency=%.3e' % (eV1, BeVortex.trans, SiO2Vortex.trans, SiVortex.absrp, net))
304 return net
307def Assemble_Collimator(eV1, xHe=1, xAl=0,xKapton=0, WD=6.0, xsw=0):
308 # from sample surface to detector
309 # xsw=0/1/-1, 6cm, 3cm, no collimator
310 # He_path depends on the collliator.
311 if xHe==1:
312 He_path=WD
313 else:
314 He_path=0.
315 air_path = WD-He_path
316 kapton_inside=0 # kapton inside collimator. no collimator, no kapton_inside
317 if xsw==0: # WD60mm collimator
318 kapton_inside=1 # collimator has thin Kapton on the second aperture
319 if xHe==1:
320 air_path = 1.51 # 1.51cm between second aperture and Be of Vortex
321 else:
322 air_path=WD
323 He_path = WD-air_path # He is between sample surface to the second aperture of xrm collimator
324 if xsw==1: # WD30mm colllimator
325 kapton_inside=1 # collimator has thin Kapton on the second aperture
326 if xHe==1: # modified 11/18/2010
327 He_path = 1.088 # from sample surface to Kapton cover/collimator
328 air_path = WD - He_path # 1.912cm between second aperture and Be of Vortex
329 xKapton = xKapton+1 # one Kapton film used to fill He from sample surface to collimator
330 else:
331 air_path = WD0; He_path=0
332 air=Material('N1.56O0.48C0.03Ar0.01Kr0.000001Xe0.0000009', 0.0013, air_path)
333 kapton=Material('C22H10O4N2', 1.42, 0.000762) # 0.3mil thick
334 HeGas=Material('He', 0.00009, He_path)
335 AlFoil=Material('Al', 2.72, 0.00381) # 1.5mil thick
336 kaptonCollimator=Material('C22H10O4N2', 1.42, 0.000762) # 0.3mil thick
337 #
338 air.getLa(eV1, 1) # 1 means number of layers here.
339 kapton.getLa(eV1, xKapton) #number of Kapton layers before collimator
340 HeGas.getLa(eV1, xHe) # number of He gas layers, default=0
341 AlFoil.getLa(eV1, xAl) # number of Al foil, default=0
342 kaptonCollimator.getLa(eV1, kapton_inside) # without collimator, no addition kapton inside
343 #
344 net=air.trans*HeGas.trans
345 net=net*kapton.trans*AlFoil.trans*kaptonCollimator.trans
346 if print2screen:
347 print('%.3f eV: air.trans=%.3e, HeGas.trans=%.3e, kapton.trans=%.3e, AlFoil.trans=%.3e, kaptonCollimator.trans=%.3e, net=%.3e' % \
348 (eV1, air.trans, HeGas.trans, kapton.trans, AlFoil.trans,kaptonCollimator.trans, net))
349 #print('%.3f eV: air.la=%.3e, HeGas.la=%.3e, kapton.la=%.3e' % (eV1, air.la, HeGas.la, kapton.la))
350 return net
353def Assemble_Detector(eV1, xHe=1, xAl=0,xKapton=0, WD=6.0, xsw=0):
354 det_efficiency=Assemble_QuadVortex(eV1)
355 trans2det=Assemble_Collimator(eV1, xHe, xAl,xKapton, WD, xsw)
356 net=det_efficiency*trans2det
357 return net
361'''
362----------------------------------------------------------------------------------------------------------
363Combine detector, transmission, sample self-absorption, and fluorescing element distribution/concentration.
364----------------------------------------------------------------------------------------------------------
365'''
366# this function is for XSW/XRM.
367def cal_NetYield2(eV0, Atoms, xHe=0, xAl=0, xKapton=0, WD=6.0, xsw=0, WriteFile='Y' , xsect_resonant=0.0, sample=''):
368 # incident energy, list of elements, experimental conditions
369 # this one tries Ka, Kb, Lg, Lb, La, Lb
370 angle0=45.; textOut=''
371 if xsw!=0 and WriteFile=='Y': print( 'XSW measurements')
372 if sample=='':
373 Include_SelfAbsorption='No'
374 else:
375 Include_SelfAbsorption='Yes'
376 angle0=sample.angle/math.pi*180.
377 textOut=sample.txt # substrate concentration
378 NetYield=[]; NetTrans=[]; Net=[]; out2=''; net=0.0
379 text0=''
380 edges =['K', 'K', 'L1', 'L1', 'L2', 'L2', 'L2', 'L3', 'L3', 'L3']
381 Fluo_lines =['Ka', 'Kb', 'Lb', 'Lg', 'Ln', 'Lb', 'Lg', 'Ll', 'La', 'Lb']
382 outputfile='Elemental_Sensitivity.txt'
383 if WriteFile=='Y':
384 fo=open(outputfile, 'w')
385 desc=' '
386 if xHe==0: desc=' not '
387 if xsw==-1:
388 out1='# 13IDC XRM/XSW using QuadVortex, incident x-ray energy at '+str(eV0)+' eV at '+str(angle0)+' degrees \n'
389 out1+='# Helium path'+desc+'used, '+str(xAl)+' Al attenuators, '+str(xKapton+1)+' Kapton attenuators, '\
390 +str(WD)+' cm working distance. \n'
391 else:
392 out1='# 13IDC XRM/XSW using QuadVortex + collimator, incident x-ray energy at '+str(eV0)+' eV at '+str(angle0)+' degrees \n'
393 out1+='# Helium path'+desc+'used, '+str(xAl)+' Al attenuators, '+str(xKapton)+' Kapton attenuators, '\
394 +str(WD)+' cm working distance. \n'
395 print( out1)
396 fo.write(out1)
397 if sample!='':
398 for stuff in Atoms:
399 text1='%6.3e %s/cm^3' % (stuff.Conc*sample.Nat1, stuff.AtSym)
400 text0=text0+' '+text1
401 textOut=textOut+' '+text0 # substrate concentration + other concentrations
402 out1='# '+text0+'\n'
403 if print2screen:
404 print( out1)
405 fo.write(out1)
406 out1='%s\t%s\t%s \t%s \t%s \t%s \t%s\n' % ('atom', 'emit', 'emit_energy', 'yield', 'transmission', 'net_sensitivity', 'sensitivity*concentration')
407 if print2screen:
408 print( out1)
409 fo.write(out1)
410 for (ii, atom) in enumerate(Atoms):
411 # MN replace:
412 atnum=f1f2.AtSym2AtNum(atom.AtSym)
413 # MN replace:
414 temp=f1f2.AtNum2f1f2Xsect(atnum, eV0)
415 xsect=temp[2] # photo-electric cross-section for fluorescence yield, temp[4] is total for attenuation
416 con=atom.Conc
417 for (nn, edge) in enumerate(edges):
418 emit=Fluo_lines[nn]
419 fy, emit_eV, emit_prob = fluo_yield(atom.AtSym, edge, emit, eV0)
420 print(emit)
421 if fy==0.0 or emit_prob==0:
422 continue # try next item if FY=0
423 else:
424 if xsect_resonant!=0.0: # use input value near edge
425 xsect=xsect_resonant # for cross-section near absoprtion edge
426 #print(xsect,fy,emit_prob)
427 net_yield=xsect*fy*emit_prob # net_yield --> cross-section, yield, emission_probability
428 # net transmission --> transmission from surface through detector
429## ------------------ [self-absorption] ------------------------------
430 trans_SelfAbsorp = 1.0 # for self-absorption.
431 if Include_SelfAbsorption=='Yes':
432 sample.getPenetration(eV0) # incident x-ray attenuation
433 if fy*emit_eV*emit_prob==0: # any one of three is zero
434 break # skip
435 sample.getLa(emit_eV)
436 # account for incident x-ray attenuation, emitted x-ray attenuation
437 trans_SelfAbsorp = sample.scale # for elements that are not part of substrate
438 if (atom in sample.ElemListFY):
439 if atom.tag=='substrate1':
440 trans_SelfAbsorp = sample.scale1 # for elements that are part of top substrate
441 if atom.tag=='substrate2':
442 trans_SelfAbsorp = sample.scale2 # for elements that are part of bottom substrate
443## ----------------------------------------------------------------------------
444 net_trans = Assemble_Detector(emit_eV, xHe, xAl,xKapton, WD, xsw)
445 net = net_yield*net_trans*trans_SelfAbsorp # elemental sensitivity
446 inten = net*con # sensitivity * concentration
447 if WriteFile=='Y':
448 out1='%s\t%s\t%6.1f \t%.3e \t%.3e \t%.3e \t%.3e\n' % (atom.AtSym+'_'+edge, emit, emit_eV, net_yield, net_trans, net, inten)
449 fo.write(out1)
450 if print2screen:
451 print('%s %s %6.1f net_yield=%.3e net_trans=%.3e net=%.3e\t' % (atom.AtSym+'_'+edge, emit, emit_eV, net_yield, net_trans, net))
452 #print(out1+' %s, depth-dependent factor= %6.4f' % (atom.tag, trans_SelfAbsorp))
453 if emit=='Kb' and fy!=0: # if above K edge, don't bother trying L edges
454 break
455 return textOut
458#def sim_spectra(eV0, Atoms, Conc, xHe=0, xAl=0, xKapton=0, WD=6.0, xsw=0, sample=''):
459def sim_spectra(eV0, Atoms, xHe=0, xAl=0, xKapton=0, WD=6.0, xsw=0, sample=''):
460 # sample=sample matrix with object attribues to add self-absorption effect
461 # Atoms is a list with elements that have attributes AtSym, Conc, tag
462 if xsw==-1: xKapton=xKapton-1 # no collimator
463 Include_SelfAbsorption='Yes'
464 Print2Screen='No'
465 if sample=='': Include_SelfAbsorption='No'
466 if xsw!=0: print('XSW measurements')
467 xx=[]; yy=[]; tag=[]; intensity_max=-10.0; LoLimit=1e-10
468 angle0=''; text1=''
469 if sample!='': # sample matrix option is used
470 angle0=str(sample.angle*180./math.pi)
471 angle0=angle0[:5]
472 for (ii, item) in enumerate(sample.ElemListFY): # add matrix elements to the lists
473 Atoms.append(item)
474 outputfile='simSpectrum_table.txt'
475 fo=open(outputfile, 'w')
476 out1='#incident x-ray at '+str(eV0)+' eV and '+angle0+' Deg.\n'
477 if Print2Screen=='Yes': print(out1)
478 fo.write(out1)
479 out1='#Emission\tenergy(eV)\tintensity \n'
480 if Print2Screen=='Yes': print(out1)
481 fo.write(out1)
482 out2='#'
483 for (ix,atom) in enumerate(Atoms):
484 # MN replace:
485 atnum=f1f2.AtSym2AtNum(atom.AtSym)
486 # con=Conc[ix]
487 con=atom.Conc
488 out2=out2+atom.AtSym+'['+str(con)+'] '
489 for edge in ['K', 'L1', 'L2', 'L3']:
490 # MN replace:
491 temp=f1f2.AtNum2f1f2Xsect(atnum, eV0)
492 xsect=temp[2] # photoelectric crosssection for each element at incident x-ray, fluorescence yield
493 # MN replace:
494 temp=elam.use_ElamFY(atom.AtSym, edge)
495 edge_eV=float(temp[2]) # absorption edge
496 if eV0>edge_eV:
497 fy=float(temp[3])
498 for EmitLine in temp[4]:
499 EmitName=EmitLine[1]
500 emit_eV=float(EmitLine[2])
501 emit_prob=float(EmitLine[3])
502 if emit_eV<fluo_emit_min: continue # ignore fluorescence below fluo_emit_min (global variable)
503 name = atom.AtSym + '_'+ EmitName
504 # net transmission --> transmission from surface through detector
505## ------------------ [self-absorption] ------------------------------
506 trans_SelfAbsorp = 1.0 # for self-absorption.
507 if Include_SelfAbsorption=='Yes':
508 sample.getPenetration(eV0) # incident x-ray attenuation
509 eV0str=str(eV0)
510 text1=' absorption_length1(%seV)= %2.2e%s \
511 absorption_length2(%seV)= %2.2e%s' % (eV0str, sample.la1*1.e4, 'microns', eV0str, sample.la2*1.e4, 'microns')
512 # text1 is added to sample.txt later
513 sample.getLa(emit_eV) # emitted fluorescence attenuation
514 trans_SelfAbsorp = sample.scale # for elements that are not part of substrate
515 if (atom in sample.ElemListFY):
516 if atom.tag=='substrate1':
517 trans_SelfAbsorp = sample.scale1 # for elements that are part of top substrate
518 if atom.tag=='substrate2':
519 trans_SelfAbsorp = sample.scale2 # for elements that are part of bottom substrate
520## ----------------------------------------------------------------------------
521 trans = Assemble_Detector(emit_eV, xHe, xAl,xKapton, WD, xsw)
522 intensity = con * fy * emit_prob * xsect * trans * trans_SelfAbsorp
523 if intensity<LoLimit: continue # skip weak emission, arbtraray limit =1e-10.
524 if intensity>intensity_max: intensity_max=intensity
525 xx.append(emit_eV); yy.append(intensity); tag.append(name)
526 for ix in range(len(yy)):
527 yy[ix]=yy[ix]/intensity_max*100.00 # makes the strongest line to 100.0
528 out1='%s\t%f\t%f \n' % (tag[ix], xx[ix], yy[ix])
529 if Print2Screen=='Yes': print(out1)
530 fo.write(out1)
531 if Print2Screen=='Yes': print(out2)
532 fo.write(out2)
533 fo.close()
534 out1=sim_GaussPeaks(xx, yy, det_res, eV0) # det_res: detector resoultion for Gaussian width (global variable)
535 if Include_SelfAbsorption=='Yes':
536 sample.txt=sample.txt+text1 # sample.txt is combined to output of cal_NetYield
537 text=cal_NetYield2(eV0, Atoms, xHe, xAl, xKapton, WD, xsw, sample=sample) # calculate net yield with weight-averaged emission
538 print(out2)
539 return text # cal_NetYield2 output is str with number densities of elements
542def sim_GaussPeaks(xx, yy, width, eV0): #xx, yy: lists, width: a peak width, eV0: incident energy
543 xline=[]; yline=[]; dX=10.0; minX=fluo_emit_min #10eV steps, lowest-->fluo_emit_min
544 amp=100 # arbitrary multiplier to shift up spectrum
545 NumOfSteps = int((eV0-minX)/dX)
546 NumOfPeaks = len(xx)
547 for ix in range(NumOfSteps+1):
548 xline.append(0); yline.append(0)
549 for (iy,peak) in enumerate(xx):
550 X0=float(peak)
551 Y0=float(yy[iy])
552 for ix in range(NumOfSteps+1):
553 energy = minX + ix*dX
554 inten = Y0*math.exp(-((energy-X0)/width)**2)*amp
555 xline[ix]=energy
556 yline[ix]+=inten
557 outputfile='simSpectrum_plot.txt'
558 fo=open(outputfile, 'w')
559# DetectorLimit=1e5 # upperlimit for total counts to 1e5 (1e5 CPS)
560 LoLimit=0.001 # low limit for each channel
561# total=0.0
562 factor=1.0
563 for ix in range(NumOfSteps+1):
564 if yline[ix]<LoLimit: yline[ix]=LoLimit
565# total+=yline[ix] # add counts
566# factor=total/DetectorLimit
567 for ix in range(NumOfSteps+1):
568 yline[ix]=yline[ix]/factor
569 out1=str(xline[ix])+'\t'+str(yline[ix])+'\n'
570 fo.write(out1)
571 fo.close()
572 #return xline, yline
575class input_param:
576 def __init__(self, eV0=14000,
577 Atoms=[],
578 xHe=0.0,
579 xAl=0,
580 xKap=0,
581 WD=6.0,
582 xsw=0):
583 if Atoms==[]:
584 atom1=ElemFY()
585 Atoms.append(atom1)
586 self.eV0=eV0
587 #incident x-ray energy in eV
588 list1=[]; list2=[]
589 for item in Atoms:
590 list1.append(item.AtSym)
591 list2.append(item.Conc)
592 self.Atoms=list1
593 #elements list
594 self.Conc=list2
595 #relative concentrations list
596 self.xHe=xHe
597 #He gas path before collimator?
598 self.xAl=xAl
599 #number of Al foils as attenuator
600 self.xKap=xKap
601 #number of Kapton foils as attenuator
602 self.WD=WD
603 #working distance in cm
604 self.xsw=xsw
605 #x-ray standing wave setup? WD=3.0 for xsw=1
609# ----------------------------------------------------------------
613if __name__=='__main__':
614 testing = 0
615 if testing:
616 #atom0='Fe'
617 #emission0='Ka'
618 #eV0=10000.
619 #edge0='K'
620 #eV1=6400.
621 #mat0=Material('SiO2', 2.2, 0.001)
622 #mat0.getLa(8000)
623 #print(mat0.la, mat0.trans)
624 #print(Assemble_QuadVortex(eV1))
625 #print(Assemble_Collimator(eV1, xHe=1, xAl=0,xKapton=0, WD=6.0, xsw=0))
626 # MN replace:
627 matrix=SampleMatrix2('CaCO3', f1f2.nominal_density('CaCO3'), 0.001,'Fe2O3', f1f2.nominal_density('Fe2O3'), 0.001, 45., 'all')
628 eV0=7500.
629 Atoms=[]
630 atom=ElemFY('La', 10e-6); Atoms.append(atom)
631 atom=ElemFY('Ce', 10e-6); Atoms.append(atom)
632 atom=ElemFY('Nd', 10e-6); Atoms.append(atom)
633 for (ii, item) in enumerate(matrix.ElemListFY):
634 pass
635 sim_spectra(eV0, Atoms, sample=matrix) # xKapton=-1 to remove WD60, -2 for WD30
636 else:
637 # March 2012 for Nov2011 beamtime
638 eV0 = 10000.
639 print2screen=0
640 Atoms = []
641 atom=ElemFY('Al', 10e-6); Atoms.append(atom)
642 atom=ElemFY('Si', 10e-6); Atoms.append(atom)
643 atom=ElemFY('Ca', 10e-6); Atoms.append(atom)
644 atom=ElemFY('Cr', 10e-6); Atoms.append(atom)
645 atom=ElemFY('Mn', 10e-6); Atoms.append(atom)
646 atom=ElemFY('Fe', 10e-6); Atoms.append(atom)
647 atom=ElemFY('Ni', 10e-6); Atoms.append(atom)
648 #sim_spectra(eV0, Atoms, xHe=1, xAl=0, xKapton=0, WD=3.0, xsw=1, sample='')
649 #Assemble_QuadVortex(1486.56)
650 #Assemble_Collimator(1486.56, 1, 0, 0, 3., 1)
651 #
652 air_path = 1.088
653 He_path = 1.912
654 air=Material('N1.56O0.48C0.03Ar0.01Kr0.000001Xe0.0000009', 0.0013, air_path)
655 kapton=Material('C22H10O4N2', 1.42, 0.000762) # 0.3mil thick
656 HeGas=Material('He', 0.00009, He_path)
657 AlFoil=Material('Al', 2.72, 0.00381) # 1.5mil thick
658 #
659 from math import *
660 emitE = [
661 1486.4, 1557.0, 1739.6, 1837.0, 3691.1, 4013.1, 5411.6,
662 5947.0, 5896.5, 6492.0, 6400.8, 7059.6, 7474.4, 8266.6
663 ]
664 for eV1 in emitE:
665 air.getLa(eV1, 1) # 1 means number of layers here.
666 print( eV1, air.la, exp(-1./air.la))
669'''
670def AtSym2FY(AtSym, Shell): #returns fluorescence yield using atomic symbol and edge(K,L1,L2,L3,M)
671 out=elam.use_ElamFY(AtSym, Shell)
672 #AtSym, edge, edge-energy, FY, [ [transition, emission, energy, probability],[]...]
673 FY=out[3]
674 return FY
676class Element:
677 def __init__(self, AtNum, AtWt, f1, f2, Xsection):
678 self.AtNum=AtNum
679 self.AtWt=AtWt #atomic weight g/mol
680 self.f1=f1 #real part of scattering factor
681 self.f2=f2 #imaginary part of scattering factor
682 self.Xsection=Xsection #atomic crossection Barns/atom
683 # f1, f2, Xsection depends on incident x-ray energy
686def setElement(AtSym, shell, eV0):
687 AtNum=AtSym2AtNum(AtSym)
688 AtWt=AtSym2AtWt(AtSym)
689 f1f2=AtNum2f1f2Xsect(AtNum, eV0)
690 AtSym=Element(AtNum, AtWt, f1f2[0], f1f2[1], f1f2[2])
691 return AtSym
694class input_param:
695 def __init__(self, eV0=14000,
696 Atoms=[],
697 xHe=0.0,
698 xAl=0,
699 xKap=0,
700 WD=6.0,
701 xsw=0):
702 if Atoms==[]:
703 atom1=ElemFY()
704 Atoms.append(atom1)
705 self.eV0=eV0
706 #incident x-ray energy in eV
707 list1=[]; list2=[]
708 for item in Atoms:
709 list1.append(item.AtSym)
710 list2.append(item.Conc)
711 self.Atoms=list1
712 #elements list
713 self.Conc=list2
714 #relative concentrations list
715 self.xHe=xHe
716 #He gas path before collimator?
717 self.xAl=xAl
718 #number of Al foils as attenuator
719 self.xKap=xKap
720 #number of Kapton foils as attenuator
721 self.WD=WD
722 #working distance in cm
723 self.xsw=xsw
724 #x-ray standing wave setup? WD=3.0 for xsw=1
725'''