/
IndirectTransmissionMonitor.py
240 lines (184 loc) · 9.13 KB
/
IndirectTransmissionMonitor.py
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
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
# Mantid Repository : https://github.com/mantidproject/mantid
#
# Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
# NScD Oak Ridge National Laboratory, European Spallation Source,
# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
# SPDX - License - Identifier: GPL - 3.0 +
#pylint: disable=no-init
from mantid.simpleapi import *
from mantid.api import *
from mantid.kernel import *
from mantid import logger
import numpy
class IndirectTransmissionMonitor(PythonAlgorithm):
_sample_ws_in = None
_can_ws_in = None
_out_ws = None
def category(self):
return "Workflow\\Inelastic;Inelastic\\Indirect"
def summary(self):
return "Calculates the sample transmission using the raw data files of the sample and its background or container."
def PyInit(self):
self.declareProperty(WorkspaceProperty('SampleWorkspace', '', direction=Direction.Input),
doc='Sample workspace')
self.declareProperty(WorkspaceProperty('CanWorkspace', '', direction=Direction.Input),
doc='Background/can workspace')
self.declareProperty(WorkspaceProperty('OutputWorkspace', '', direction=Direction.Output),
doc='Output workspace group')
def PyExec(self):
setup_prog = Progress(self, start=0.0, end=0.05, nreports=5)
setup_prog.report('Setting up algorithm')
self._setup()
ws_basename = str(self._sample_ws_in)
trans_prog = Progress(self, start=0.05, end=0.4, nreports=2)
trans_prog.report('Transforming monitor for Sample')
self._trans_mon(ws_basename, 'Sam', self._sample_ws_in)
trans_prog.report('Transforming monitor for Container')
self._trans_mon(ws_basename, 'Can', self._can_ws_in)
workflow_prog = Progress(self, start=0.4, end=1.0, nreports=4)
# Generate workspace names
sam_ws = ws_basename + '_Sam'
can_ws = ws_basename + '_Can'
trans_ws = ws_basename + '_Trans'
# Divide sample and can workspaces
workflow_prog.report('Dividing Sample by Container')
Divide(LHSWorkspace=sam_ws, RHSWorkspace=can_ws, OutputWorkspace=trans_ws)
trans = numpy.average(mtd[trans_ws].readY(0))
logger.information('Average Transmission: ' + str(trans))
workflow_prog.report('Adding Sample logs')
AddSampleLog(Workspace=trans_ws, LogName='can_workspace', LogType='String', LogText=self._can_ws_in)
# Group workspaces
workflow_prog.report('Creating output GroupWorkspace')
group = sam_ws + ',' + can_ws + ',' + trans_ws
GroupWorkspaces(InputWorkspaces=group, OutputWorkspace=self._out_ws)
self.setProperty('OutputWorkspace', self._out_ws)
workflow_prog.report('Algorithm complete')
def _setup(self):
"""
Get properties.
"""
self._sample_ws_in = self.getPropertyValue("SampleWorkspace")
self._can_ws_in = self.getPropertyValue("CanWorkspace")
self._out_ws = self.getPropertyValue('OutputWorkspace')
def _get_spectra_index(self, input_ws):
"""
Gets the index of the two monitors and first detector for the current instrument configurtion.
Assumes monitors are named monitor1 and monitor2
"""
workspace = mtd[input_ws]
instrument = workspace.getInstrument()
# Get workspace index of first detector
detector_1_idx = 2
try:
# First try to get first detector for current analyser bank
analyser = instrument.getStringParameter('analyser')[0]
detector_1_idx = instrument.getComponentByName(analyser)[0].getID() - 1
logger.information('Got index of first detector for analyser %s: %d' % (analyser, detector_1_idx))
except IndexError:
# If that fails just get the first spectrum which is a detector
spectrumInfo = workspace.spectrumInfo()
for spec_idx in range(workspace.getNumberHistograms()):
if not spectrumInfo.isMonitor(spec_idx):
detector_1_idx = spec_idx
logger.information('Got index of first detector in workspace: %d' % (detector_1_idx))
break
# Get workspace index of monitor(s)
monitor_1_idx = 0
monitor_2_idx = None
if instrument.hasParameter('Workflow.Monitor1-SpectrumNumber'):
# First try to get monitors based on workflow parameters in IPF
monitor_1_idx = int(instrument.getNumberParameter('Workflow.Monitor1-SpectrumNumber')[0])
if instrument.hasParameter('Workflow.Monitor2-SpectrumNumber'):
monitor_2_idx = int(instrument.getNumberParameter('Workflow.Monitor2-SpectrumNumber')[0])
logger.information('Got index of monitors: %d, %s' % (monitor_1_idx, str(monitor_2_idx)))
else:
# If that fails just use some default values (correct ~60% of the time)
monitor_2_idx = 1
logger.warning('Could not determine index of monitors, using default values.')
return monitor_1_idx, monitor_2_idx, detector_1_idx
def _get_detector_workspace_index(self, workspace, detector_id):
"""
Returns the workspace index for a given detector ID in a workspace.
@param workspace Workspace to find detector in
@param detector_id Detector ID to search for
"""
for spec_idx in range(0, mtd[workspace].getNumberHistograms()):
if mtd[workspace].getDetector(spec_idx).getID() == detector_id:
return spec_idx
return None
def _unwrap_mon(self, input_ws):
out_ws = '_unwrap_mon_out'
# Unwrap monitor - inWS contains M1,M2,S1 - outWS contains unwrapped Mon
# Unwrap s1>2 to L of S2 (M2) ie 38.76 Ouput is in wavelength
_, join = UnwrapMonitor(InputWorkspace=input_ws, OutputWorkspace=out_ws, LRef='37.86')
# Fill bad (dip) in spectrum
RemoveBins(InputWorkspace=out_ws,
OutputWorkspace=out_ws,
Xmin=join-0.001,
Xmax=join+0.001,
Interpolation="Linear")
FFTSmooth(InputWorkspace=out_ws,
OutputWorkspace=out_ws,
WorkspaceIndex=0,
IgnoreXBins=True)
DeleteWorkspace(input_ws)
return out_ws
def _trans_mon(self, ws_basename, file_type, input_ws):
monitor_1_idx, monitor_2_idx, detector_1_idx = self._get_spectra_index(input_ws)
CropWorkspace(InputWorkspace=input_ws,
OutputWorkspace='__m1',
StartWorkspaceIndex=monitor_1_idx,
EndWorkspaceIndex=monitor_1_idx)
if monitor_2_idx is not None:
CropWorkspace(InputWorkspace=input_ws,
OutputWorkspace='__m2',
StartWorkspaceIndex=monitor_2_idx,
EndWorkspaceIndex=monitor_2_idx)
CropWorkspace(InputWorkspace=input_ws,
OutputWorkspace='__det',
StartWorkspaceIndex=detector_1_idx,
EndWorkspaceIndex=detector_1_idx)
# Check for single or multiple time regimes
mon_tcb_start = mtd['__m1'].readX(0)[0]
spec_tcb_start = mtd['__det'].readX(0)[0]
DeleteWorkspace('__det')
mon_ws = '__Mon'
if spec_tcb_start == mon_tcb_start:
mon_ws = self._unwrap_mon('__m1') # unwrap the monitor spectrum and convert to wavelength
RenameWorkspace(InputWorkspace=mon_ws,
OutputWorkspace='__Mon1')
else:
ConvertUnits(InputWorkspace='__m1',
OutputWorkspace='__Mon1',
Target="Wavelength")
mon_ws = ws_basename + '_' + file_type
if monitor_2_idx is not None:
ConvertUnits(InputWorkspace='__m2',
OutputWorkspace='__Mon2',
Target="Wavelength")
DeleteWorkspace('__m2')
x_in = mtd['__Mon1'].readX(0)
xmin1 = mtd['__Mon1'].readX(0)[0]
xmax1 = mtd['__Mon1'].readX(0)[len(x_in) - 1]
x_in = mtd['__Mon2'].readX(0)
xmin2 = mtd['__Mon2'].readX(0)[0]
xmax2 = mtd['__Mon2'].readX(0)[len(x_in) - 1]
wmin = max(xmin1, xmin2)
wmax = min(xmax1, xmax2)
CropWorkspace(InputWorkspace='__Mon1',
OutputWorkspace='__Mon1',
XMin=wmin,
XMax=wmax)
RebinToWorkspace(WorkspaceToRebin='__Mon2',
WorkspaceToMatch='__Mon1',
OutputWorkspace='__Mon2')
Divide(LHSWorkspace='__Mon2',
RHSWorkspace='__Mon1',
OutputWorkspace=mon_ws)
DeleteWorkspace('__Mon1')
DeleteWorkspace('__Mon2')
else:
RenameWorkspace(InputWorkspace='__Mon1',
OutputWorkspace=mon_ws)
# Register algorithm with Mantid
AlgorithmFactory.subscribe(IndirectTransmissionMonitor)