Newer
Older
from MantidFramework import *
from mantidsimple import *
import math
# Keep track of whether we have numpy since readX() may
# return a numpy array or a list depending on whether it
# is available.
try:
import numpy
HAS_NUMPY = True
except:
class EQSANSTransmission(PythonAlgorithm):
"""
Calculate transmission for EQSANS
Transmission counts is the signal going through a pinhole in the middle of the
beam stop (after hitting the sample). The number of counts is equal to the neutron
flux times the transmission N_trans = T*Phi. The transmission T is obtained by normalizing
the N_trans by the monitor or accelerator current. T = N_trans/current.
Dividing the pixel counts by T we account for transmission.
"""
# Size of the box around the beam center
NX_TRANS_PIX = 10
NY_TRANS_PIX = 10
# To define transmission peak
transmission_peak_to_bg_ratio = 1000
# To define transmission peak mask
min_transmission_peak_to_bg_ratio = 5
def category(self):
return "SANS"
def name(self):
return "EQSANSTransmission"
def PyInit(self):
# Input workspace
self.declareWorkspaceProperty("InputWorkspace", "", Direction.Input)
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
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
# Output workspace to put the transmission histo into
self.declareProperty("OutputWorkspace", "")
# X position of the beam center
self.declareProperty("XCenter", 96.0)
# Y position of the beam center
self.declareProperty("YCenter", 128.0)
# Transmission will be normalized to 1 if True
self.declareProperty("NormalizeToUnity", False)
def PyExec(self):
input_ws = self.getProperty("InputWorkspace")
output_ws = self.getProperty("OutputWorkspace")
xcenter = int(math.floor(self.getProperty("XCenter")))
ycenter = int(math.floor(self.getProperty("YCenter")))
normalize = self.getProperty("NormalizeToUnity")
nx_min = xcenter - self.NX_TRANS_PIX
nx_max = xcenter + self.NX_TRANS_PIX
ny_min = ycenter - self.NY_TRANS_PIX
ny_max = ycenter + self.NY_TRANS_PIX
# Sum up all TOF bins
Integration(input_ws, input_ws+'_int')
# Find pixel with highest and lowest signal
counts_2D = numpy.zeros([2*self.NY_TRANS_PIX+1, 2*self.NX_TRANS_PIX+1])
# Smallest count in a given pixel in the region around the center
signal_min = None
# Highest count in a given pixel in the region around the center
signal_max = None
# Position of the pixel with the highest count
xmax = xcenter
ymax = ycenter
for ix in range(nx_min, nx_max+1):
for iy in range(ny_min, ny_max+1):
i_pixel = 256*ix+iy
signal = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel)[0]
if signal_min is None or signal < signal_min:
signal_min = signal
if signal_max is None or signal > signal_max:
signal_max = signal
xmax = ix
ymax = iy
self.log().information("Transmission search [%d, %d] [%d, %d]" % (nx_min, nx_max, ny_min, ny_max))
self.log().information("Max transmission counts at [%d, %d] = %g" % (xmax, ymax,signal_max))
self.log().information("Min transmission counts = %g" % (signal_min))
# Find transmission peak boundaries
# The peak is defined by the area around the pixel with the highest count
# that is above the background level.
xpeakmax = nx_max
for i in range(xmax+1, nx_max+1):
i_pixel = 256*i+ymax
signal = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel)[0]
if signal*self.transmission_peak_to_bg_ratio < signal_max:
xpeakmax = i-1
break
i_pixel_low = 256*(i-1)+ymax
signal_low = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel_low)[0]
if signal > signal_low:
xpeakmax = i-1
break
xpeakmin = nx_min
for i in range(xmax-1, nx_min+1, -1):
i_pixel = 256*i+ymax
signal = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel)[0]
if signal*self.transmission_peak_to_bg_ratio < signal_max:
xpeakmin = i+1
break
i_pixel_high = 256*(i+1)+ymax
signal_high = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel_high)[0]
if signal > signal_high:
xpeakmin = i+1
break
ypeakmax = ny_max
for i in range(ymax+1, ny_max+1):
i_pixel = 256*xmax+i
signal = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel)[0]
if signal*self.transmission_peak_to_bg_ratio < signal_max:
ypeakmax = i-1
break
i_pixel_low = 256*xmax+(i-1)
signal_low = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel_low)[0]
if signal > signal_low:
ypeakmax = i-1
break
ypeakmin = ny_min
for i in range(ymax-1, ny_min+1, -1):
i_pixel = 256*xmax+i
signal = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel)[0]
if signal*self.transmission_peak_to_bg_ratio < signal_max:
ypeakmin = i+1
break
i_pixel_high = 256*xmax+(i-1)
signal_high = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel_high)[0]
if signal > signal_high:
ypeakmin = i+1
break
self.log().information("Peak range in X [%d, %d]" % (xpeakmin, xpeakmax))
self.log().information("Peak range in Y [%d, %d]" % (ypeakmin, ypeakmax))
# Mask monitor peak around the beam center
# Question for JK: why do we do this?
mask = []
for ix in range(xpeakmin, xpeakmax+1):
for iy in range(ypeakmin, ypeakmax+1):
i_pixel = 256*ix+iy
mask.append(i_pixel)
for iy in range(ypeakmin, ypeakmax+1):
for ix in range(xpeakmin, xcenter):
i_pixel = 256*ix+iy
signal = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel)[0]
i_pixel_high = 256*(ix+1)+iy
signal_high = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel_high)[0]
if signal > signal_high and signal*self.min_transmission_peak_to_bg_ratio<signal_max:
if i_pixel in mask:
mask.remove(i_pixel)
else:
break
for iy in range(ypeakmin, ypeakmax+1):
for ix in range(xpeakmax, xcenter, -1):
i_pixel = 256*ix+iy
signal = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel)[0]
i_pixel_high = 256*(ix-1)+iy
signal_high = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel_high)[0]
if signal > signal_high and signal*self.min_transmission_peak_to_bg_ratio<signal_max:
if i_pixel in mask:
mask.remove(i_pixel)
else:
break
for iy in range(xpeakmin, xpeakmax+1):
for ix in range(ypeakmin, ycenter):
i_pixel = 256*ix+iy
signal = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel)[0]
i_pixel_high = 256*ix+iy+1
signal_high = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel_high)[0]
if signal > signal_high and signal*self.min_transmission_peak_to_bg_ratio<signal_max:
if i_pixel in mask:
mask.remove(i_pixel)
else:
break
for iy in range(xpeakmin, xpeakmax+1):
for ix in range(ypeakmax, ycenter, -1):
i_pixel = 256*ix+iy
signal = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel)[0]
i_pixel_high = 256*ix+iy-1
signal_high = mantid.getMatrixWorkspace(input_ws+'_int').readY(i_pixel_high)[0]
if signal > signal_high and signal*self.min_transmission_peak_to_bg_ratio<signal_max:
if i_pixel in mask:
mask.remove(i_pixel)
else:
break
# Sum transmission peak counts as a function of TOF
dataX = mantid.getMatrixWorkspace(input_ws).readX(0)
if HAS_NUMPY:
dataX = dataX.tolist()
nTOF = len(dataX)
dataY = (nTOF-1)*[0]
dataE = (nTOF-1)*[0]
transmission_counts = 0
total_pixels = 0
# We will normalize the transmission by the accelerator current
proton_charge = mantid.getMatrixWorkspace(input_ws).getRun()["proton_charge"].getStatistics().mean
duration = mantid.getMatrixWorkspace(input_ws).getRun()["proton_charge"].getStatistics().duration
frequency = mantid.getMatrixWorkspace(input_ws).getRun()["frequency"].getStatistics().mean
acc_current = 1.0e-12 * proton_charge * duration * frequency
for itof in range(nTOF-1):
for ix in range(xpeakmin, xpeakmax+1):
for iy in range(ypeakmin, ypeakmax+1):
i_pixel = 256*ix+iy
signal = mantid.getMatrixWorkspace(input_ws).readY(i_pixel)[itof]
error = mantid.getMatrixWorkspace(input_ws).readE(i_pixel)[itof]
if i_pixel in mask:
total_pixels += 1
dataY[itof] += signal
dataE[itof] += error*error
transmission_counts += signal
dataE[itof] = math.sqrt(dataE[itof])/acc_current
dataY[itof] /= acc_current
if normalize:
self.log().information("Normalizing the translation to average 1")
factor = transmission_counts/acc_current/(nTOF-1)
for itof in range(nTOF-1):
dataY[itof] /= factor
dataE[itof] /= factor
unitX = mantid.getMatrixWorkspace(input_ws).getAxis(0).getUnit().name()
CreateWorkspace(output_ws, dataX, dataY, dataE, NSpec=1, UnitX=unitX)
total_pixels /= (nTOF-1)
self.log().information("Total transmission counts (%d) = %g" % (total_pixels, transmission_counts))
mtd.registerPyAlgorithm(EQSANSTransmission())