cena_energy_resolution_kazama¶
Analyze Kazama-san’s energy resolution data
The program input is the numeric data from Figs 16-18 in Kazama-2006, derived using “g3data”.
'''Analyze Kazama-san's energy resolution data
The program input is the numeric data from Figs 16-18 in Kazama-2006,
derived using "g3data".
'''
import os
import sys
import logging
logging.basicConfig()
import datetime
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import scipy.optimize
matrix_g3data_debug = True
def matrix_g3data(visualize=False):
''' Graced data matrix in Kazama-2006.
:returns: Three branches, 25, 750, and 3300 eV, is returned, i.e.,
(matrix_25, matrix_750, matrix_3300).
matrix_`XX` is (N, 2) numpy array, where N is the data point.
Data is retrieved using ``g3data`` manually.
Then, the data are interpolated.
'''
# 25 eV three data
dat25_1 = np.fromstring('''2.64150943396 0
7.54716981132 34.1523341523
12.641509434 122.604422604
17.5471698113 516.216216216
22.4528301887 947.42014742
27.5471698113 1004.91400491
32.641509434 929.72972973
37.5471698113 743.98034398
42.4528301887 617.936117936
47.5471698113 487.469287469
52.4528301887 339.312039312
57.5471698113 184.520884521
62.641509434 124.815724816
67.5471698113 80.5896805897
72.8301886792 47.4201474201
77.5471698113 38.5749385749
82.641509434 18.6732186732
87.5471698113 9.82800982801
92.4528301887 9.82800982801
97.5471698113 7.61670761671
102.641509434 3.1941031941
107.547169811 0.982800982801''', sep=' ').reshape((22, 2))
dat25_2 = np.fromstring('''2.64150943396 0
7.54716981132 34.353993788
12.641509434 122.693430995
17.5471698113 513.985443419
22.4528301887 945.080895647
27.5471698113 1004.67340411
32.641509434 929.376477678
37.5471698113 743.518612953
42.4528301887 615.154605721
47.5471698113 484.575124009
52.4528301887 338.520699087
57.5471698113 183.616893051
62.4528301887 128.225858792
67.5471698113 79.4645588985
72.641509434 46.1823744843
77.5471698113 37.2286866626
82.641509434 15.0030133049
87.5471698113 8.26062769459
92.4528301887 8.15214871819
97.5471698113 5.82819526216
102.641509434 1.29293959483
107.547169811 1.18446061842
''', sep=' ').reshape((22, 2))
dat25_3 = np.fromstring('''2.53343823761 0
7.44295830055 34.2485360337
12.4468922109 121.656248821
17.5452399685 522.347223472
22.4547600315 948.501194689
27.4586939418 1006.01960858
32.4626278521 932.910716204
37.4665617624 744.672672533
42.5649095201 613.998333532
47.474429583 486.646769694
52.4783634933 337.154113477
57.4822974036 185.44743512
62.4862313139 126.729686652
67.4901652242 81.2960710252
72.4940991345 46.9325660999
77.4036191975 38.0311867635
82.5019669552 18.0579547731
87.5059008655 8.04869339016
92.5098347758 8.00253163822
97.5137686861 7.95636988628
102.517702596 2.37515278379
107.521636507 1.22197996174''', sep=' ').reshape((22, 2))
dat25 = np.zeros_like(dat25_1)
dat25[:, 0] = np.linspace(2.5, 107.5, num=22)
dat25[:, 1] = np.array([dat25_1[:, 1], dat25_2[:, 1], dat25_3[:, 1]]).mean(0) / 10000
dat750_1 = np.fromstring('''250 0
354.722202882 201.10701107
449.63621806 322.878228782
549.51959897 498.15498155
649.751096568 747.232472325
759.407853512 994.464944649
858.377428114 976.014760148
951.811947365 784.132841328
1050.39859361 684.501845018
1148.83729026 553.505535055
1252.06259138 437.269372694
1350.62312887 332.103321033
1453.83972708 214.022140221
1552.72227251 177.121771218
1651.47427418 112.546125461
1750.37422544 79.3357933579
1849.29158254 49.815498155
1948.27856297 35.0553505535
2051.96511871 16.6051660517
2151.02172248 16.6051660517
2250.0522175 11.0701107011
2353.79969366 5.53505535055
2448.1393163 5.53505535055
2551.89549537 1.84501845018
2655.67778319 3.69003690037
''', sep=' ').reshape((25, 2))
dat750_2 = np.fromstring('''250.196695515 0
349.331235248 198.339483395
450.826121164 328.413284133
549.960660897 500.922509225
649.095200629 760.147601476
750.590086546 993.542435424
849.724626279 976.014760148
951.219512195 778.597785978
1050.35405193 678.966789668
1149.48859166 547.04797048
1250.98347758 431.734317343
1350.11801731 330.258302583
1449.25255704 214.022140221
1550.74744296 177.121771218
1649.88198269 113.468634686
1749.01652242 78.4132841328
1850.51140834 49.815498155
1949.64594807 34.1328413284
2051.14083399 15.6826568266
2150.27537372 15.6826568266
2249.40991345 9.22509225092
2350.90479937 5.53505535055
2450.0393391 4.61254612546
2551.53422502 0
2653.02911094 1.84501845018
''', sep=' ').reshape((25, 2))
dat750_3 = np.fromstring('''
250.505050505 0.2457002457
349.494949495 194.656019656
450.841750842 325.491400491
549.831649832 498.71007371
651.178451178 759.459459459
750.168350168 993.488943489
851.515151515 975.982800983
950.505050505 781.572481572
1049.49494949 682.063882064
1150.84175084 546.621621622
1249.83164983 435.135135135
1351.17845118 330.098280098
1450.16835017 214.926289926
1551.51515152 178.071253071
1650.50505051 113.574938575
1749.49494949 79.484029484
1850.84175084 50
1949.83164983 35.257985258
2051.17845118 16.8304668305
2150.16835017 16.8304668305
2251.51515152 10.3808353808
2352.86195286 6.69533169533
2444.78114478 5.77395577396
2553.1986532 1.16707616708
2649.83164983 3.9312039312
''', sep=' ').reshape((25, 2))
dat750 = np.zeros_like(dat750_1)
dat750[:, 0] = np.linspace(250, 2650, num=25)
dat750[:, 1] = np.array([dat750_1[:, 1], dat750_2[:, 1], dat750_3[:, 1]]).mean(0) / 10000
dat3300_1 = np.fromstring('''
1376.03305785 0
1623.96694215 18.8747731397
1871.90082645 178.584392015
2119.83471074 429.764065336
2380.16528926 543.012704174
2628.09917355 694.010889292
2876.03305785 707.078039927
3123.96694215 734.664246824
3371.90082645 650.453720508
3619.83471074 577.858439201
3867.76859504 463.157894737
4128.09917355 399.274047187
4376.03305785 294.736842105
4623.96694215 200.362976407
4871.90082645 162.613430127
5119.83471074 85.6624319419
5380.16528926 75.499092559
5628.09917355 34.8457350272
5876.03305785 26.1343012704
6123.96694215 18.8747731397
6371.90082645 5.80762250454
6619.83471074 8.71143375681
6880.16528926 10.1633393829
7128.09917355 4.3557168784
7376.03305785 7.25952813067
7623.96694215 2.90381125227
7871.90082645 4.3557168784
''', sep=' ').reshape((27, 2))
dat3300_2 = np.fromstring('''
1376.03305785 0
1623.96694215 18.8747731397
1871.90082645 178.584392015
2119.83471074 429.764065336
2380.16528926 543.012704174
2628.09917355 694.010889292
2876.03305785 707.078039927
3123.96694215 734.664246824
3371.90082645 650.453720508
3619.83471074 577.858439201
3867.76859504 463.157894737
4128.09917355 399.274047187
4376.03305785 294.736842105
4623.96694215 200.362976407
4871.90082645 162.613430127
5119.83471074 85.6624319419
5380.16528926 75.499092559
5628.09917355 34.8457350272
5876.03305785 26.1343012704
6123.96694215 18.8747731397
6371.90082645 5.80762250454
6619.83471074 8.71143375681
6880.16528926 10.1633393829
7128.09917355 4.3557168784
7376.03305785 7.25952813067
7623.96694215 2.90381125227
7871.90082645 4.3557168784
''', sep=' ').reshape((27, 2))
dat3300_3 = np.fromstring('''
1376.43932684 0.726744186047
1624.44641275 17.4418604651
1872.45349867 180.23255814
2126.66076174 431.686046512
2374.66784765 535.610465116
2622.67493357 694.76744186
2876.88219663 707.122093023
3124.88928255 735.465116279
3372.89636847 650.436046512
3627.10363153 575.581395349
3875.11071745 462.936046512
4123.11780337 401.162790698
4377.32506643 294.331395349
4625.33215235 199.854651163
4873.33923826 162.063953488
5127.54650133 84.3023255814
5375.55358725 75.5813953488
5623.56067316 34.8837209302
5877.76793623 25.4360465116
6125.77502214 18.8953488372
6373.78210806 5.08720930233
6627.98937112 7.99418604651
6875.99645704 9.4476744186
7124.00354296 3.63372093023
7378.21080602 6.54069767442
7626.21789194 2.90697674419
7874.22497786 3.63372093023
''', sep=' ').reshape((27, 2))
dat3300 = np.zeros_like(dat3300_1)
dat3300[:, 0] = np.linspace(1375, 7875, num=27)
dat3300[:, 1] = np.array([dat3300_1[:, 1], dat3300_2[:, 1], dat3300_3[:, 1]]).mean(0) / 10000
if visualize:
fig = plt.figure()
enestep = np.linspace(0, 3, 300)
ax = fig.add_subplot(221)
ax.plot(dat25_1[:, 0], dat25_1[:, 1])
ax.plot(dat25_2[:, 0], dat25_2[:, 1])
ax.plot(dat25_3[:, 0], dat25_3[:, 1])
ax.plot(dat25[:, 0], dat25[:, 1] * 1e4, lw=3)
lim0 = 12.5
lim1 = 62.5
ax.plot(enestep * 25., energy_response_triangle(enestep, limit0=lim0/25., limit1=lim1/25.) * dat25[:, 1].max() * 1.1e4)
ax.text(80, 800, '%g--%g eV' % (lim0, lim1))
ax.set_title('25 eV calibration')
ax = fig.add_subplot(222)
ax.plot(dat750_1[:, 0], dat750_1[:, 1])
ax.plot(dat750_2[:, 0], dat750_2[:, 1])
ax.plot(dat750_3[:, 0], dat750_3[:, 1])
ax.plot(dat750[:, 0], dat750[:, 1] * 1e4, lw=3)
lim0 = 325
lim1 = 1650
ax.plot(enestep * 750., energy_response_triangle(enestep, limit0=lim0/750., limit1=lim1/750.) * dat750[:, 1].max() * 1e4)
ax.text(2000, 800, '%g--%g eV' % (lim0, lim1))
ax.set_title('750 eV calibration')
ax = fig.add_subplot(223)
ax.plot(dat3300_1[:, 0], dat3300_1[:, 1])
ax.plot(dat3300_2[:, 0], dat3300_2[:, 1])
ax.plot(dat3300_3[:, 0], dat3300_3[:, 1])
ax.plot(dat3300[:, 0], dat3300[:, 1] * 1e4, lw=3)
peak = 3000.
lim0 = 1650.
lim1 = 4950.
ax.plot(enestep * 3300., energy_response_triangle(enestep, peak=peak/3300., limit0=lim0/3300, limit1=lim1/3300.) * dat3300[:, 1].max() * 1.1e4)
ax.text(4000, 800, '%g--(%g)--%g eV' % (lim0, peak, lim1))
ax.set_title('3300 eV calibration')
# Normalized
ax = fig.add_subplot(224)
ax.plot(dat25[:, 0] / 25.0, dat25[:, 1] / dat25[:, 1].max(), label='25eV')
ax.plot(dat750[:, 0] / 750.0, dat750[:, 1] / dat750[:, 1].max(), label='750eV')
ax.plot(dat3300[:, 0] / 3300.0, dat3300[:, 1] / dat3300[:, 1].max(), label='3300eV')
ax.plot(enestep, energy_response_triangle(enestep), ':', lw=2, label='Calrep.')
ax.grid(True)
ax.legend()
return dat25, dat750, dat3300
def step_function(x, n):
''' A step function
:param x: Scalar or np.array.
:param n: Step function parameter.
:returns: 0 if x<n whereas 1 for x>=n.
'''
x2 = np.array(x)
mask = (x2 >= n)
x2 = np.ma.zeros(x2.shape)
x2.mask = mask
return x2.filled(1) * 1.
def energy_response_triangle(e, peak=1.0, limit0=0.2, limit1=2.2):
''' Energy response function with triangle shape
According to Martin's calibration report, triangle approximation with
half max in 0.6--1.6 E0 is implemented.
This means that E<=0.2 and E>=2.2 is no response.
:param e: Energy normalized by the incoming beam energy
:type e: ``np.array``
:returns: Energy response for the given energy, *e*, normalized by the
maximum resonse at E=1. (energy_response_triangle(1) == 1)
'''
e2 = np.array(e)
e2 = np.ma.masked_less(e2, limit0)
e2 = np.ma.masked_greater(e2, limit1)
ke = (((e2 - peak) / (peak - limit0) + 1) * (1. - step_function(e2, peak))
+ ((e2 - peak) / (peak - limit1) + 1) * step_function(e2, peak))
return ke.filled(0) * 1. # * 1. is for scalar input to return as scalar
if __name__ == '__main__':
mx = matrix_g3data(visualize=True)