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

1""" 

2Fluorescence Intensity calculations: 

3 Y. Choi 

4 

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' 

15 

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""" 

45 

46 

47import math 

48import numpy 

49import sys 

50from larch.utils.physical_constants import AVOGADRO, BARN 

51from larch.xray import xray_delta_beta, chemparse, atomic_mass 

52 

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 

56 

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 

87 

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 

97 

98 

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 

105 

106 

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 

277 

278 

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. 

291 

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 

305 

306 

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 

351 

352 

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 

358 

359 

360 

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 

456 

457 

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 

540 

541 

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 

573 

574 

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 

606 

607 

608 

609# ---------------------------------------------------------------- 

610 

611 

612 

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)) 

667 

668 

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 

675 

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 

684 

685 

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 

692 

693 

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'''