forked from AmbaPant/mantid
-
Notifications
You must be signed in to change notification settings - Fork 1
/
VesuvioAnalysis.py
315 lines (279 loc) · 17.5 KB
/
VesuvioAnalysis.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
# 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 +
'''
This file was last edited by Giovanni Romanelli on 23/05/2020.
The structure of the script is as follows:
First, the collection of functions allowing the reduction and analysis of the spectra;
Second, a list of input parameters specific of the VESUVIO experiment;
Last, the reduction and analysis procedure, iterative in TOF and finally in y-space for hydrogen.
The script has been tested to be run in Mantid Workbench in a Windows operative system.
PLEASE, DO NOT MODIFY THE "TECHNICAL SECTION" UNLESS YOU ARE AN
EXPERT VESUVIO INSTRUMENT SCIENTIST.
'''
##########################################################
#### TECHNICAL SECTION - NOT FOR USERS
##########################################################
from __future__ import (absolute_import, division, print_function, unicode_literals)
from mantid.kernel import *
from mantid.api import *
from mantid.simpleapi import *
from mantid.kernel import StringListValidator, IntListValidator, FloatBoundedValidator, Direction
from Inelastic.vesuvio.analysisHelpers import *
################################################################################################
##################
##################
##################
################## HAVE F
##################
##################
##################
##################
################################################################################################
##########################################################
#### USER SECTION - FOR USERS
##########################################################
'''
The user section is composed of an initialisation section, an iterative analysis/reduction section
of the spectra in the time-of-flight domain, and a final section where the analysis of the corrected
hydrogen neutron Compton profile is possible in the Y-space domain.
The fit procedure in the time-of-flight domain is based on the scipy.minimize.optimize() tool,
used with the SLSQP minimizer, that can handle both boundaries and constraints for fitting parameters.
The Y-space analysis is, at present, performed on a single spectrum, being the result of
the sum of all the corrected spectra, subsequently symmetrised and unit-area normalised.
The Y-space fit is performed using the Mantid minimiser and average Mantid resolution function, using
a Gauss-Hermite expansion including H0 and H4 at present, while H3 (proportional to final-state effects)
is not needed as a result of the symmetrisation.
'''
class VesuvioAnalysis(PythonAlgorithm):
def category(self):
return "Inelastic\\Indirect\\Vesuvio"
def seeAlso(self):
return ["VesuvioCalculateGammaBackground", "ConvertToYSpace","VesuvioThickness",
"VesuvioCalculateMS","Integration","VesuvioResolution","Fit",]
def PyInit(self):
self.declareProperty(
"AnalysisMode",
"LoadReduceAnalyse",
doc="In the first case, all the algorithm is run. In the second case, the data are not re-loaded, and only"
" the TOF and y-scaling bits are run. In the third case, only the y-scaling final analysis is run.",
validator=StringListValidator(
[
"LoadReduceAnalyse",
"ReduceAnalyse",
"Analyse"]))
self.declareProperty(
FileProperty(
"IPFile",
"ip2018.par",
action=FileAction.Load,
direction=Direction.Input,
extensions=["par"]),
doc="The instrument parameter file")
self.declareProperty("NumberOfIterations",2, doc="Number of time the reduction is reiterated.",
validator=IntListValidator([0,1,2,3,4]))
self.declareProperty("OutputName", "polyethylene",doc="The base name for the outputs." )
self.declareProperty("Runs", "38898-38906", doc="List of Vesuvio run numbers (e.g. 20934-20937, 30924)")
self.declareProperty(IntArrayProperty("Spectra",[135,182]), doc="Range of spectra to be analysed (first, last).")
self.declareProperty(FloatArrayProperty("TOFRangeVector", [110.,1.5,460.]),
doc="In micro seconds (lower bound, binning, upper bound).")
self.declareProperty(
"TransmissionGuess",
0.9174,
doc="A number from 0 to 1 to represent the experimental transmission value of the sample for epithermal"
" neutrons. This value is used for the multiple scattering corrections. If 1, the multiple scattering correction is not run.",
validator=FloatBoundedValidator(
0,
1) )
self.declareProperty("MultipleScatteringOrder", 2, doc="Order of multiple scattering events in MC simultation.",
validator=IntListValidator([1,2,3,4]))
self.declareProperty("MonteCarloEvents", 1.e6, doc="Number of events for MC multiple scattering simulation.")
self.declareProperty(ITableWorkspaceProperty("ComptonProfile",
"",
direction=Direction.Input),doc="Table for Compton profiles")
self.declareProperty(IntArrayProperty("ConstraintsProfileNumbers", [0,1]))
self.declareProperty(
"ConstraintsProfileScatteringCrossSection",
"2.*82.03/5.551",
doc="The ratio of the first to second intensities, each equal to atom stoichiometry times bound scattering"
"cross section. Simple arithmetic can be included but the result may be rounded.")
self.declareProperty("ConstraintsProfileState", "eq", validator=StringListValidator(["eq","ineq"]))
self.declareProperty(IntArrayProperty("SpectraToBeMasked", []))#173,174,181
self.declareProperty("SubtractResonancesFunction", "", doc="Function for resonance subtraction. Empty means no subtraction.")
self.declareProperty("YSpaceFitFunctionTies", "",doc="The TOF spectra are subtracted by all the fitted profiles"
" about the first element specified in the elements string. Then such spectra are converted to the Y space"
" of the first element (using the ConvertToYSPace algorithm). The spectra are summed together and"
" symmetrised. A fit on the resulting spectrum is performed using a Gauss Hermte function up to the sixth"
"order.")
def validateInputs(self):
tableCols = [
"symbol",
"mass(a.u.)",
"Intensity lower limit",
"Intensity value",
"Intensity upper limit",
"Width lower limit",
"Width value",
"Width upper limit",
"Centre lower limit",
"Centre value",
"Centre upper limit"]
issues = dict()
table = self.getProperty("ComptonProfile").value
if(not table):
issues["ComptonProfile"] = "An elements table should be provided."
elif(table.columnCount()!= len(tableCols) or sorted(cleanNames(tableCols))!=sorted(cleanNames(table.getColumnNames())) ):
issues["ComptonProfile"] = "The table should be of the form: "
for name in tableCols:
issues["ComptonProfile"] += name + ", "
TOF = self.getProperty("TOFRangeVector").value
if len(TOF) != 3:
issues["TOFRangeVector"] = "TOFRangeVector should have length 3 (lower, binning, upper)."
constraints = self.getProperty("ConstraintsProfileNumbers").value
if len(constraints)!=2:
issues["ConstraintsProfileNumbers"] = "ConstraintsProfileNumbers should only contain 2 numbers."
#check aritmatic is safe
cross_section = self.getProperty("ConstraintsProfileScatteringCrossSection").value
for ch in cross_section:
if ch not in ["+","-","*","/",".","(",")"] and not ch.isdigit():
issues["ConstraintsProfileScatteringCrossSection"]= "Must be a valid mathmatical expression. "+ch
spectra = self.getProperty("Spectra").value
if len(spectra) != 2:
issues["Spectra"] = "Spectra should be of the form [first, last]"
run_string = self.getProperty("Runs").value
for ch in run_string:
if ch not in ["-",","] and not ch.isdigit():
issues["Runs"]= "Runs are list, can use - for a range or commas to seperate runs. "
masks = self.getProperty("SpectraToBeMasked").value
spec = [j for j in range(spectra[0], spectra[1]+1)]
for mask in masks:
if mask not in spec:
issues["SpectraToBeMasked"] = "Masked spectra is not in the loaded spectra."
return issues
def PyExec(self):
IPFile = self.getProperty("IPFile").value
g_log = logger(self.log())
analysisMode = self.getProperty("AnalysisMode").value
# This is the number of iterations for the reduction analysis in time-of-flight.
number_of_iterations = self.getProperty("NumberOfIterations").value
# Parameters of the measurement
ws_name=self.getProperty("OutputName").value
runs = self.getProperty("Runs").value
first_spectrum, last_spectrum = self.getProperty("Spectra").value
tof_range = self.getProperty("TOFRangeVector").value
# Parameters for the multiple-scattering correction, including the shape of the sample.
transmission_guess = self.getProperty("TransmissionGuess").value
multiple_scattering_order = self.getProperty("MultipleScatteringOrder").value
number_of_events = self.getProperty("MonteCarloEvents").value
# parameters of the neutron Compton profiles to be fitted.
elements = generate_elements(self.getProperty("ComptonProfile").value)
# constraint on the intensities of element peaks
#provide LHS element, RHS element, mult. factor, flag
# if flag=True inequality; if flag = False equality
constraints_profile_num = self.getProperty("ConstraintsProfileNumbers").value
# check this is valid in the validate inputs
cross_section = evaluate(self.getProperty("ConstraintsProfileScatteringCrossSection").value)
state = self.getProperty("ConstraintsProfileState").value
C1 = constraint( constraints_profile_num[0], constraints_profile_num[1], cross_section ,state)
constraints = [C1]
# spectra to be masked
spectra_to_be_masked = self.getProperty("SpectraToBeMasked").value
subtract_resonances = True
resonance_function = self.getProperty("SubtractResonancesFunction").value
if not resonance_function:
subtract_resonances = False
fit_hydrogen_in_Y_space = True# If True, corrected time-of-flight spectra containing H only are transformed to Y-space and fitted.
y_fit_ties = self.getProperty("YSpaceFitFunctionTies").value
if not y_fit_ties:
fit_hydrogen_in_Y_space = False
####### END OF USER INPUT
#
# Start of the reduction and analysis procedure
#
if "Load" in analysisMode:
spectrum_list=str(first_spectrum)+'-'+str(last_spectrum)
LoadVesuvio(
Filename=runs,
SpectrumList=spectrum_list,
Mode="SingleDifference",
SumSpectra=False,
InstrumentParFile=IPFile,
OutputWorkspace=ws_name)
# chose limits such that there is at list one non-nan bin among all the spectra between -30 and 30 \AA-1
Rebin(InputWorkspace=ws_name,Params=tof_range,OutputWorkspace=ws_name)
if "Reduce" in analysisMode:
vertical_width, horizontal_width, thickness = 0.1, 0.1, 0.001 # expressed in meters
create_slab_geometry(ws_name,vertical_width, horizontal_width, thickness)
masses, par, bounds, constraints = prepare_fit_arguments(elements, constraints)
fit_arguments = [bounds, constraints]
# Iterative analysis and correction of time-of-flight spectra.
for iteration in range(number_of_iterations):
if iteration == 0:
ws_to_be_fitted = CloneWorkspace(InputWorkspace = ws_name, OutputWorkspace = ws_name+"_cor")
ws_to_be_fitted = mtd[ws_name+"_cor"]
MaskDetectors(Workspace=ws_to_be_fitted,SpectraList=spectra_to_be_masked)
# Fit and plot where the spectra for the current iteration
spectra, widths, intensities, centres = block_fit_ncp(
par, first_spectrum,last_spectrum, masses, ws_to_be_fitted, fit_arguments, g_log,IPFile, g_log)
# Calculate mean widths and intensities
mean_widths, mean_intensity_ratios = calculate_mean_widths_and_intensities(
masses, widths, intensities, spectra, g_log) # at present is not multiplying for 0,9
if (number_of_iterations - iteration -1 > 0):
# evaluate gamma background correction ---------- This creates a
# background workspace with name : str(ws_name)+"_gamma_background"
sample_properties = calculate_sample_properties(masses, mean_widths, mean_intensity_ratios, "GammaBackground", g_log)
correct_for_gamma_background(ws_name, first_spectrum,last_spectrum, sample_properties, g_log)#
Scale(InputWorkspace = str(ws_name)+"_gamma_background", OutputWorkspace = str(ws_name)+"_gamma_background",
Factor=0.9, Operation = "Multiply")
# evaluate multiple scattering correction --------- This creates a
# background workspace with name : str(ws_name)+"_MulScattering"
if transmission_guess < 1. :
sample_properties = calculate_sample_properties(
masses, mean_widths, mean_intensity_ratios, "MultipleScattering", g_log)
correct_for_multiple_scattering(ws_name, first_spectrum,last_spectrum, sample_properties, transmission_guess,
multiple_scattering_order, number_of_events,g_log, masses,mean_intensity_ratios)
# Create corrected workspace
Minus(LHSWorkspace= ws_name, RHSWorkspace = str(ws_name)+"_gamma_background",
OutputWorkspace = ws_name+"_cor")
if transmission_guess < 1. :
Minus(LHSWorkspace= ws_name+"_cor", RHSWorkspace = str(ws_name)+"_MulScattering",
OutputWorkspace = ws_name+"_cor")
if subtract_resonances :
Minus(LHSWorkspace=ws_name+"_cor",RHSWorkspace=ws_name+"_cor_fit",
OutputWorkspace=ws_name+"_cor_residuals")
ws = CloneWorkspace(ws_name+"_cor_residuals")
for index in range( ws.getNumberHistograms() ):
Fit(Function=resonance_function, InputWorkspace=ws_name+"_cor_residuals", WorkspaceIndex=index,
MaxIterations=10000,
Output=ws_name+"_cor_residuals", OutputCompositeMembers=True, StartX=110., EndX=460.)
fit_ws=mtd[ws_name+"_cor_residuals_Workspace"]
for bin in range( ws.blocksize() ) :
ws.dataY(index)[bin]=fit_ws.readY(1)[bin]
ws.dataE(index)[bin]=0.
RenameWorkspace(InputWorkspace="ws", OutputWorkspace=ws_name+"_fitted_resonances")
Minus(LHSWorkspace=ws_name+"_cor",RHSWorkspace=ws_name+"_fitted_resonances",
OutputWorkspace=ws_name+"_cor" )
Minus(LHSWorkspace=ws_name+"_cor_residuals",RHSWorkspace=ws_name+"_fitted_resonances",
OutputWorkspace=ws_name+"_cor_residuals" )
else:
if fit_hydrogen_in_Y_space:
hydrogen_ws = subtract_other_masses(ws_name+"_cor", widths, intensities, centres, spectra, masses,IPFile,g_log)
RenameWorkspace(hydrogen_ws, ws_name+'_H')
SumSpectra(InputWorkspace=ws_name+'_H', OutputWorkspace=ws_name+'_H_sum')
calculate_mantid_resolutions(ws_name, masses[0])
# Fit of the summed and symmetrised hydrogen neutron Compton profile in its Y space using MANTID.
if fit_hydrogen_in_Y_space:
#calculate_mantid_resolutions(ws_name, masses[0])
#max_Y
_ = convert_to_y_space_and_symmetrise(ws_name+"_H",masses[0])
# IT WOULD BE GOOD TO HAVE THE TIES DEFINED IN THE USER SECTION!!!
constraints = " sigma1=3.0,c4=0.0, c6=0.0,A=0.08, B0=0.00, ties = {}".format(y_fit_ties)
correct_for_offsets = True
y_range = (-20., 20.)
final_fit(ws_name+'_H_JoY_sym', constraints, y_range,correct_for_offsets, masses, g_log)
elif "Analyse" in analysisMode:
g_log.notice("Did not compute analysis. The YSpaceFitFunctionTies must be stated.")
AlgorithmFactory.subscribe(VesuvioAnalysis)