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
|
from __future__ import division, print_function
import numpy as np
from bumps.names import Parameter, pmath, FitProblem, cosd, sind
from .peaks import Peaks, Gaussian, Background
def read_data():
# data= Z1.T
X = np.linspace(0.4840, 0.5080, 13)
Y = np.linspace(-0.5180, -0.4720, 23)
X, Y = np.meshgrid(X, Y)
A = np.genfromtxt('XY_mesh2.txt', unpack=True)
Z1 = A[26:39]
data = Z1.T
err = np.sqrt(data)
#err= A[39:54]
return X, Y, data, err
def build_problem():
M = Peaks([Gaussian(name="G1-"),
Gaussian(name="G2-"),
#Gaussian(name="G3-"),
#Gaussian(name="G4-"),
Background()],
*read_data())
background = np.min(M.data)
background += np.sqrt(background)
signal = np.sum(M.data) - M.data.size*background
M.parts[-1].C.value = background
peak1 = M.parts[0]
peak1.xc.range(0.45, 0.55)
peak1.yc.range(-0.55, -0.4)
peak1.xc.value = 0.500
peak1.yc.value = -0.485
if 0:
# Peak centers are independent
for peak in M.parts[1:-1]:
peak.xc.range(0.45, 0.55)
peak.yc.range(-0.55, -0.4)
M.parts[1].xc.value = 0.495
M.parts[1].yc.value = -0.495
else:
# Peak centers lie on a line
theta = Parameter(45, name="theta")
theta.range(0, 90)
for i, peak in enumerate(M.parts[1:-1]):
delta = Parameter(0.0045, name="delta-%d"%(i+1))
delta.range(0.0, 0.015)
peak.xc = peak1.xc + delta*cosd(theta)
peak.yc = peak1.yc + delta*sind(theta)
# Initial values
cx, cy = 0.4996-0.4957, -0.4849+0.4917
theta.value = np.degrees(np.arctan2(cy, cx))
delta.value = np.sqrt(cx**2 + cy**2)
# Initial values
for peak in M.parts[:-1]:
peak.A.value = signal/(len(M.parts)-1) # Equal size peaks
dx, dy = 0.4997-0.4903, -0.4969+0.4851
dxm, dym = 0.4951-0.4960, -0.4941+0.4879
peak1.s1.value = np.sqrt(dx**2 + dy**2)/2.35/2
peak1.s2.value = np.sqrt(dxm**2 + dym**2)/2.35/2
peak1.theta.value = np.degrees(np.arctan2(dy, dx))
# Peak intensity varies
for peak in M.parts[:-1]:
peak.A.range(0.1*signal, 1.1*signal)
peak1.s1.range(0.002, 0.02)
peak1.s2.range(0.001, 0.02)
peak1.theta.range(-90, -0)
if 1:
# Peak shape is the same across all peaks
for peak in M.parts[1:-1]:
peak.s1 = peak1.s1
peak.s2 = peak1.s2
peak.theta = peak1.theta
else:
for peak in M.parts[1:-1]:
peak.s1.range(*peak1.s1.bounds.limits)
peak.s2.range(*peak1.s2.bounds.limits)
peak.theta.range(*peak1.theta.bounds.limits)
if 1:
M.parts[-1].C.pmp(100.0)
if 1:
for peak in M.parts[:-1]:
peak.s1.value = 0.006
peak.s2.value = 0.002
peak.theta.value = -60.0
peak.A.value = signal/2
if 0:
print("shape", peak1.s1.value, peak1.s2.value, peak1.theta.value)
print("centers theta,delta", theta.value, delta.value)
print("centers", (peak1.xc.value, peak1.yc.value),
(M.parts[1].xc.value, M.parts[1].yc.value))
return FitProblem(M)
problem = build_problem()
# Note: if you want to run the model file as a script without invoking bumps
# (which can be handy for debugging data loaders, etc.) then add the script
# statements to a *__name__ == "__main__"* code block at the end of the file.
# For example:
#
# if __name__ == "__main__":
# print("data", read_data())
#
# Since this model file imports from the *.peaks* module, you
# also need to set up the package context to resolve the dot in
# the relative import statement. This can be done by adding the
# following two lines before the .peaks import at the top of the file:
#
# from bumps.names import relative_import
# __package__ = relative_import(__file__)
|
