-
Notifications
You must be signed in to change notification settings - Fork 122
/
PEARLTransfit.py
366 lines (342 loc) · 18.1 KB
/
PEARLTransfit.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
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
# Mantid Repository : https://github.com/mantidproject/mantid
#
# Copyright © 2020 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 +
# -----------------------------------------------------------------------
# Transfit v2 - translated from Fortran77 to Python/Mantid
# Original Author: Christopher J Ridley
# Last updated: 2nd October 2020
#
# Program fits a Voigt function to PEARL transmission data, using the doppler broadening
# of the gaussian component to determine the sample temperature
#
#
# References:
# 'Temperature measurement in a Paris-Edinburgh cell by neutron resonance spectroscopy' - Journal Of Applied Physics 98,
# 064905 (2005)
# 'Remote determination of sample temperature by neutron resonance spectroscopy' - Nuclear Instruments and Methods in
# Physics Research A 547 (2005) 601-615
# -----------------------------------------------------------------------
from mantid.kernel import Direction, StringListValidator
from mantid.api import PythonAlgorithm, MultipleFileProperty, AlgorithmFactory, mtd
from mantid.simpleapi import (LoadRaw, DeleteWorkspace, CropWorkspace, NormaliseByCurrent, ExtractSingleSpectrum,
ConvertUnits, CreateWorkspace, Fit, Plus, Multiply, RenameWorkspace,
CreateSingleValuedWorkspace)
from os import path
from scipy import constants
import numpy as np
class PEARLTransfit(PythonAlgorithm):
# Resonance constants broadly consistent with those reported from
# 'Neutron cross sections Vol 1, Resonance Parameters'
# S.F. Mughabghab and D.I. Garber
# June 1973, BNL report 325
# Mass in atomic units, En in eV, temperatures in K, gamma factors in eV
ResParamsDict = {
# Hf01
"Hf01_Mass": 177.0,
"Hf01_En": 1.098,
"Hf01_TD": 252.0,
"Hf01_TwogG": 0.00192,
"Hf01_Gg": 0.0662,
"Hf01_startE": 0.6,
"Hf01_Ediv": 0.01,
"Hf01_endE": 1.7,
# Hf02
"Hf02_Mass": 177.0,
"Hf02_En": 2.388,
"Hf02_TD": 252.0,
"Hf02_TwogG": 0.009,
"Hf02_Gg": 0.0608,
"Hf02_startE": 2.0,
"Hf02_Ediv": 0.01,
"Hf02_endE": 2.7,
# Ta10
"Ta10_Mass": 181.0,
"Ta10_En": 10.44,
"Ta10_TD": 240.0,
"Ta10_TwogG": 0.00335,
"Ta10_Gg": 0.0069,
"Ta10_startE": 9.6,
"Ta10_Ediv": 0.01,
"Ta10_endE": 11.4,
# Irp6
"Irp6_Mass": 191.0,
"Irp6_En": 0.6528,
"Irp6_TD": 420.0,
"Irp6_TwogG": 0.000547,
"Irp6_Gg": 0.072,
"Irp6_startE": 0.1,
"Irp6_Ediv": 0.01,
"Irp6_endE": 0.9,
# Iro5
"Iro5_Mass": 191.0,
"Iro5_En": 5.36,
"Iro5_TD": 420.0,
"Iro5_TwogG": 0.006,
"Iro5_Gg": 0.082,
"Iro5_startE": 4.9,
"Iro5_Ediv": 0.01,
"Iro5_endE": 6.3,
# Iro9
"Iro9_Mass": 191.0,
"Iro9_En": 9.3,
"Iro9_TD": 420.0,
"Iro9_TwogG": 0.0031,
"Iro9_Gg": 0.082,
"Iro9_startE": 8.7,
"Iro9_Ediv": 0.01,
"Iro9_endE": 9.85
}
# Physical constants
k = constants.k
e = constants.e
Pi = np.pi
def version(self):
return 1
def name(self):
return 'PEARLTransfit'
def category(self):
return 'Diffraction\\Fitting'
def summary(self):
return "Reads high-energy neutron resonances from the downstream monitor data on the PEARL instrument," \
"then fits a Voigt function to them to determine the sample temperature. A calibration must be run for" \
" each sample pressure. Can be used on a single file, or multiple files, in which case workspaces" \
" are summed and the average taken."
def PyInit(self):
self.declareProperty(MultipleFileProperty('Files', extensions=[".raw", ".s0x"]),
doc='Files of calibration runs (numors). Must be detector scans.')
self.declareProperty(name='FoilType',
defaultValue='Hf01',
validator=StringListValidator(['Hf01', 'Hf02', 'Ta10', 'Irp6', 'Iro5', 'Iro9']),
direction=Direction.Input,
doc="Type of foil included with the sample")
self.declareProperty(name='Ediv',
defaultValue='0.0025',
direction=Direction.Input,
doc="Energy step in eV")
self.declareProperty(name='ReferenceTemp',
defaultValue='290',
direction=Direction.Input,
doc="Enter reference temperature in K")
self.declareProperty(name='Calibration',
defaultValue=False,
direction=Direction.Input,
doc="Calibration flag, default is False in which case temperature measured")
self.declareProperty(name='Debug',
defaultValue=False,
direction=Direction.Input,
doc="True/False - provides more verbose output of procedure for debugging purposes")
def validateFileInputs(self, filesList):
# MultipleFileProperty returns a list of list(s) of files, or a list containing a single (string) file
# Condense into one list of file(s) in either case
returnList = []
for files in filesList:
if isinstance(files, str):
returnList.append(files)
else:
returnList.extend(files)
return returnList
def PyExec(self):
# ----------------------------------------------------------
# Imports resonance parameters from ResParam dictionary depending on the selected foil type.
# ----------------------------------------------------------
files = self.getProperty("Files").value
files = self.validateFileInputs(files)
foilType = self.getProperty("FoilType").value
divE = float(self.getProperty("Ediv").value)
isCalib = self.getProperty("Calibration").value
mass = self.ResParamsDict[foilType + '_Mass']
TD = self.ResParamsDict[foilType + '_TD']
# Energy parameters are in eV
energy = self.ResParamsDict[foilType + '_En']
TwogG = self.ResParamsDict[foilType+'_TwogG']
Gg = self.ResParamsDict[foilType+'_Gg']
startE = self.ResParamsDict[foilType+'_startE']
endE = self.ResParamsDict[foilType+'_endE']
refTemp = float(self.getProperty("ReferenceTemp").value)
isDebug = self.getProperty("Debug").value
if not isCalib:
if 'S_fit_Parameters' not in mtd:
self.log().warning("No calibration files found. Please run this algorithm will 'Calibration' ticked to "
"generate the calibration workspace.")
return
self.monitorTransfit(files, foilType, divE)
file = files[0]
discard, fileName = path.split(file)
fnNoExt = path.splitext(fileName)[0]
# Define the gaussian width at the reference temperature
# Factor of 1e3 converting to meV for Mantid
width_300 = 1000.0 * np.sqrt(4 * energy * self.k * refTemp * mass / (self.e * (1 + mass) ** 2))
# ----------------------------------------------------------
# Perform fits based on whether isCalib is flagged (calibration or measurement)
# ----------------------------------------------------------
if isCalib:
lorentzFWHM = 1000.0 * (0.5 * TwogG + Gg)
# Take peak position starting guess from tabulated value
peakPosGuess = energy * 1000
fileName_monitor = fnNoExt + '_monitor'
# For guessing initial background values, read in y-data and use starting value as value for b0
wsBgGuess = mtd[fileName_monitor]
bg0guess = 0.86 * wsBgGuess.readY(0)[0]
# New Voigt function as from Igor pro function
Fit(Function='name=PEARLTransVoigt,Position=' + str(peakPosGuess) + ',LorentzianFWHM=' + str(
lorentzFWHM) + ',GaussianFWHM=' + str(width_300) + ',Amplitude=1.6,Bg0=' + str(
bg0guess) + ',Bg1=0.0063252,Bg2=0,constraints=(1<LorentzianFWHM,' + str(width_300)
+ '<GaussianFWHM)', InputWorkspace=fileName_monitor, MaxIterations=200, Output='S_fit')
#
DeleteWorkspace('S_fit_NormalisedCovarianceMatrix')
DeleteWorkspace(fileName_monitor)
else:
S_fit = mtd['S_fit_Parameters']
lorentzFWHM = S_fit.column(1)[1]
gaussianFWHM = S_fit.column(1)[2]
# Calculates the gaussian resolution contribution, sometimes constraints are broken and the value can drop
# below that from width_300 alone, in which case the instrument contribution is set to zero.
if gaussianFWHM > width_300:
gaussianFWHM_inst = np.sqrt(gaussianFWHM**2 - width_300**2)
else:
gaussianFWHM_inst = 0
#
# New Voigt function as from Igor pro function
Fit(Function='name=PEARLTransVoigt,Position=' + str(S_fit.column(1)[0]) + ',LorentzianFWHM='
+ str(lorentzFWHM) + ',GaussianFWHM=' + str(gaussianFWHM) + ',Amplitude='
+ str(S_fit.column(1)[3]) + ',Bg0=' + str(S_fit.column(1)[4]) + ',Bg1='
+ str(S_fit.column(1)[5]) + ',Bg2=' + str(S_fit.column(1)[6]) + ',constraints=('
+ str(gaussianFWHM) + '<GaussianFWHM),ties=(LorentzianFWHM=' + str(lorentzFWHM) + ')',
InputWorkspace=fnNoExt + '_monitor', MaxIterations=200, Output='T_fit')
#
DeleteWorkspace('T_fit_NormalisedCovarianceMatrix')
DeleteWorkspace(fnNoExt + '_monitor')
T_fit = mtd['T_fit_Parameters']
gaussian_FWHM_fitted = T_fit.column(1)[2]
width_T = np.sqrt(gaussian_FWHM_fitted ** 2 - gaussianFWHM_inst ** 2)
# Factor of 1e-3 converts back from meV to eV
Teff = (((width_T * 1e-3) ** 2) * self.e * ((1 + mass) ** 2)) / (4 * 1e-3 * T_fit.column(1)[0] * self.k
* mass)
Teff_low = ((((width_T - T_fit.column(2)[2]) * 1e-3) ** 2) * self.e
* ((1 + mass) ** 2))/(4 * 1e-3 * (T_fit.column(1)[0] + T_fit.column(2)[0]) * self.k * mass)
Teff_high = ((((width_T + T_fit.column(2)[2]) * 1e-3) ** 2) * self.e * ((1 + mass) ** 2)) /\
(4 * 1e-3 * (T_fit.column(1)[0] - T_fit.column(2)[0]) * self.k * mass)
errTeff = 0.5*(Teff_high-Teff_low)
# ----------------------------------------------------------
# If the temperature is too far below the Debye temperature, then the result is inaccurate. Else the
# temperature is calculated assuming free gas formulation
# ----------------------------------------------------------
if 8 * Teff < 3 * TD:
self.log().information("The effective temperature is currently too far below the Debye temperature to"
"give an accurate measure.")
Tactual = Teff
Terror = errTeff
else:
Tactual = 3 * TD / (4 * np.log((8 * Teff + 3 * TD)/(8 * Teff - 3 * TD)))
Tactual_high = 3 * TD / (4 * np.log((8 * (Teff + errTeff) + 3 * TD) / (8 * (Teff + errTeff) - 3 * TD)))
Tactual_low = 3 * TD / (4 * np.log((8 * (Teff - errTeff) + 3 * TD) / (8 * (Teff - errTeff) - 3 * TD)))
Terror = 0.5 * (np.abs(Tactual - Tactual_high) + np.abs(Tactual - Tactual_low))
# Introduce a catchment for unphysically small determined errors, setting to defualt value if below a set
# threshold
Terror_flag = 0
if Terror < 5.0:
Terror_flag = 1
Terror = 10.0
if isDebug:
self.log().information("-----------------------------")
self.log().information("Debugging....")
self.log().information("The Debye temperature is " + str(TD) + " K")
self.log().information("The effective temperature is: {:.1f}"
.format(Teff) + "+/- {:.1f}".format(errTeff) + " K")
self.log().information("Energy bin width set to " + str(1000 * divE) + " meV")
self.log().information("E range is between " + str(startE) + " and " + str(endE))
self.log().information("Gaussian width at this reference temperature is: {:.2f}".format(width_300)
+ " meV")
self.log().information("Lorentzian FWHM is fixed: {:.2f}".format(lorentzFWHM)+" meV")
self.log().information("Gaussian FWHM is fitted as: {:.2f}".format(gaussian_FWHM_fitted) + " meV")
self.log().information("Instrumental contribution is: {:.2f}".format(gaussianFWHM_inst) + " meV")
self.log().information("Temperature contribution is: {:.2f}".format(width_T) + " meV")
self.log().information("-----------------------------")
self.log().information("Sample temperature is: {:.1f}".format(Tactual) + " +/- {:.1f}".format(Terror)
+ " K")
if Terror_flag == 1:
self.log().information("(the default error, as determined error unphysically small)")
# ----------------------------------------------------------
# Define function for importing raw monitor data, summing, normalising, converting units, and cropping
# ----------------------------------------------------------
def monitorTransfit(self, files, foilType, divE):
isFirstFile = True
isSingleFile = len(files) == 1
firstFileName = ""
for file in files:
discard, fileName = path.split(file)
fnNoExt = path.splitext(fileName)[0]
if isFirstFile:
firstFileName = fnNoExt
fileName_Raw = fnNoExt + '_raw'
fileName_3 = fnNoExt + '_3'
LoadRaw(Filename=file, OutputWorkspace=fileName_Raw)
CropWorkspace(InputWorkspace=fileName_Raw, OutputWorkspace=fileName_Raw, XMin=100, XMax=19990)
NormaliseByCurrent(InputWorkspace=fileName_Raw, OutputWorkspace=fileName_Raw)
ExtractSingleSpectrum(InputWorkspace=fileName_Raw, OutputWorkspace=fileName_3, WorkspaceIndex=3)
DeleteWorkspace(fileName_Raw)
ConvertUnits(InputWorkspace=fileName_3, Target='Energy', OutputWorkspace=fileName_3)
self.TransfitRebin(fileName_3, fileName_3, foilType, divE)
if not isFirstFile:
Plus(LHSWorkspace=firstFileName + '_3', RHSWorkspace=fileName_3, OutputWorkspace=firstFileName + '_3')
DeleteWorkspace(fileName_3)
else:
isFirstFile = False
if isSingleFile:
RenameWorkspace(InputWorkspace=firstFileName + '_3', OutputWorkspace=firstFileName + '_monitor')
else:
noFiles = len(files) ** (-1)
CreateSingleValuedWorkspace(OutputWorkspace='scale', DataValue=noFiles)
Multiply(LHSWorkspace=firstFileName + '_3', RHSWorkspace='scale',
OutputWorkspace=firstFileName + '_monitor')
DeleteWorkspace('scale')
DeleteWorkspace(firstFileName + '_3')
# ----------------------------------------------------------
# Define function for improved rebinning of the data
# ----------------------------------------------------------
def TransfitRebin(self, inputWS, outputWSName, foilType, divE):
ws2D = mtd[inputWS]
# Expand the limits for rebinning to prevent potential issues at the boundaries
startE = self.ResParamsDict[foilType + '_startE']
endE = self.ResParamsDict[foilType + '_endE']
startEp = 0.99 * startE
endEp = 1.01 * endE
CropWorkspace(InputWorkspace=ws2D, OutputWorkspace=ws2D, XMin=1000 * startEp, XMax=1000 * endEp)
numPts = np.int(((endEp - startEp) / divE) + 1)
xData_out = []
xData_in = ws2D.readX(0)
yData_in = ws2D.readY(0)
# Deals with issues in Mantid where Xdata mark start of bin, not true x position
# calculates the bin width and then takes middle of bin as x value
current_bin_widths = []
xActual = []
for x in range(0, len(xData_in) - 1):
current_bin_widths.append(xData_in[x + 1] - xData_in[x])
xActual.append(xData_in[x] + 0.5 * (xData_in[x + 1] - xData_in[x]))
# Make xData with uniform binning defined by divE
for j in range(numPts):
xData_out.append(1000 * (startEp + j * divE))
yData_out = [0] * (len(xData_out))
# Normalise output ydata accordingly based on surrounding values
yNorm = [0] * (len(xData_out))
for j in range(0, len(yData_in)):
currentBin = np.int((xActual[j] - startEp * 1000) / (divE * 1000))
scale1 = 1 - ((xActual[j] - xData_out[currentBin]) / (divE * 1000))
yData_out[currentBin] += yData_in[j] * scale1
yNorm[currentBin] += scale1
if currentBin < (len(xData_out) - 1):
yData_out[currentBin + 1] += yData_in[j] * (1 - scale1)
yNorm[currentBin + 1] += 1 - scale1
# Apply the normalisation, with a catch for any potential divide by zero errors
for i in range(len(yData_out)):
if yNorm[i] != 0:
yData_out[i] = yData_out[i] / yNorm[i]
else:
print('Empty bin')
outputWS = CreateWorkspace(DataX=xData_out, DataY=yData_out, NSpec=1, UnitX='meV')
CropWorkspace(InputWorkspace=outputWS, OutputWorkspace=outputWS, XMin=1000 * startE, XMax=1000 * endE)
RenameWorkspace(InputWorkspace=outputWS, OutputWorkspace=outputWSName)
AlgorithmFactory.subscribe(PEARLTransfit)