1 | """*WIKI* |
---|
2 | |
---|
3 | |
---|
4 | *WIKI*""" |
---|
5 | |
---|
6 | from MantidFramework import * |
---|
7 | from mantidsimple import * |
---|
8 | from mantid.api import AlgorithmFactory |
---|
9 | import os |
---|
10 | |
---|
11 | all_algs = AlgorithmFactory.getRegisteredAlgorithms(True) |
---|
12 | if 'GatherWorkspaces' in all_algs: |
---|
13 | HAVE_MPI = True |
---|
14 | import boostmpi as mpi |
---|
15 | else: |
---|
16 | HAVE_MPI = False |
---|
17 | |
---|
18 | COMPRESS_TOL_TOF = .01 |
---|
19 | |
---|
20 | class SNSPowderReduction2(PythonAlgorithm): |
---|
21 | class PDConfigFile(object): |
---|
22 | class PDInfo: |
---|
23 | """Inner class for holding configuration information for a reduction.""" |
---|
24 | def __init__(self, data, has_dspace=False, has_vnoise=False): |
---|
25 | if data is None: |
---|
26 | data = [None, None, 1, 0, 0, 0., 0.] |
---|
27 | self.freq = data[0] |
---|
28 | self.wl = data[1] |
---|
29 | self.bank = int(data[2]) |
---|
30 | self.van = int(data[3]) |
---|
31 | self.can = int(data[4]) |
---|
32 | self.vnoise = 0 # default value |
---|
33 | self.has_dspace = has_dspace |
---|
34 | self.tmin = 0. # default value |
---|
35 | self.tmax = 0. # default value |
---|
36 | |
---|
37 | # calculate the remaining indices |
---|
38 | offset = 5 |
---|
39 | if has_vnoise: |
---|
40 | self.vnoise = int(data[offset]) |
---|
41 | offset += 1 |
---|
42 | |
---|
43 | if has_dspace: |
---|
44 | self.dmin = data[offset] |
---|
45 | self.dmax = data[offset+1] |
---|
46 | offset += 2 |
---|
47 | if len(data) > offset: |
---|
48 | self.tmin = data[offset] |
---|
49 | if len(data) > offset+1: |
---|
50 | self.tmax = data[offset+1] |
---|
51 | |
---|
52 | def __init__(self, filename): |
---|
53 | if len(filename.strip()) <= 0: |
---|
54 | filename = None |
---|
55 | self.filename = filename |
---|
56 | self._data = {} |
---|
57 | self.use_dspace = False |
---|
58 | self._use_vnoise = False |
---|
59 | self._focusPos = None |
---|
60 | self.iparmFile = None |
---|
61 | if self.filename is None: |
---|
62 | return |
---|
63 | handle = file(filename, 'r') |
---|
64 | lines = handle.readlines() |
---|
65 | handle.close() |
---|
66 | |
---|
67 | # create the focus positions |
---|
68 | (lines, self._focusPos) = self._generateFocusPos(lines) |
---|
69 | if len(lines) == 0: |
---|
70 | self.filename = None |
---|
71 | return |
---|
72 | |
---|
73 | # get the rest of the characterization information |
---|
74 | for line in lines: |
---|
75 | self._addData(line) |
---|
76 | |
---|
77 | def _generateFocusPos(self, lines): |
---|
78 | if not lines[0].startswith("Instrument parameter file:"): |
---|
79 | return (lines, None) |
---|
80 | |
---|
81 | result = {} |
---|
82 | |
---|
83 | # get name of parameter file |
---|
84 | temp = lines[0] |
---|
85 | temp = temp.replace("Instrument parameter file:", "") |
---|
86 | self.iparmFile = temp.strip() |
---|
87 | if len(self.iparmFile) <= 0: |
---|
88 | self.iparmFile = None |
---|
89 | lines = lines[1:] # delete this line |
---|
90 | |
---|
91 | # get the spectra into a buffer |
---|
92 | spectrainfo = [] |
---|
93 | for line in lines: |
---|
94 | if line.startswith("L1"): |
---|
95 | break |
---|
96 | spectrainfo.append(line) |
---|
97 | numSpectra = len(spectrainfo) |
---|
98 | |
---|
99 | result['PrimaryFlightPath'] = lines[numSpectra].split()[1] |
---|
100 | |
---|
101 | # delete the rest of the focus position info |
---|
102 | lines = lines[numSpectra+1:] |
---|
103 | |
---|
104 | # parse the focus positions |
---|
105 | specids = [] |
---|
106 | l2 = [] |
---|
107 | polar = [] |
---|
108 | azimuthal = [] |
---|
109 | for spec in spectrainfo: |
---|
110 | temp = spec.split() |
---|
111 | specids.append(int(temp[0])) |
---|
112 | l2.append(float(temp[1])) |
---|
113 | polar.append(float(temp[2])) |
---|
114 | azimuthal.append(0.) |
---|
115 | |
---|
116 | # assign to the correct place |
---|
117 | result['SpectrumIDs'] = specids |
---|
118 | result['L2'] = l2 |
---|
119 | result['Polar'] = polar |
---|
120 | result['Azimuthal'] = azimuthal |
---|
121 | |
---|
122 | return (lines, result) |
---|
123 | |
---|
124 | def _addData(self, line): |
---|
125 | if line.startswith('#') or len(line.strip()) <= 0: |
---|
126 | if "d_min" in line and "d_max" in line: |
---|
127 | self.use_dspace = True |
---|
128 | if "vanadium_back" in line: |
---|
129 | self._use_vnoise = True |
---|
130 | return |
---|
131 | data = line.strip().split() |
---|
132 | data = [float(i) for i in data] |
---|
133 | if data[0] not in self._data.keys(): |
---|
134 | self._data[data[0]]={} |
---|
135 | info = self.PDInfo(data, self.use_dspace, self._use_vnoise) |
---|
136 | self._data[info.freq][info.wl]=info |
---|
137 | def __getFrequency(self, request): |
---|
138 | for freq in self._data.keys(): |
---|
139 | if 100. * abs(float(freq)-request)/request < 5.: |
---|
140 | return freq |
---|
141 | raise RuntimeError("Failed to find frequency: %fHz" % request) |
---|
142 | |
---|
143 | def __getWavelength(self, frequency, request): |
---|
144 | for wavelength in self._data[frequency].keys(): |
---|
145 | if 100. * abs(wavelength-request)/request < 5.: |
---|
146 | return wavelength |
---|
147 | raise RuntimeError("Failed to find wavelength: %fAngstrom" % request) |
---|
148 | |
---|
149 | def getInfo(self, frequency, wavelength): |
---|
150 | #print "getInfo(%f, %f)" % (frequency, wavelength) |
---|
151 | if self.filename is not None: |
---|
152 | if frequency is None: |
---|
153 | raise RuntimeError("Unable to determine frequency from data") |
---|
154 | if wavelength is None: |
---|
155 | raise RuntimeError("Unable to determine wavelength from data") |
---|
156 | frequency = self.__getFrequency(float(frequency)) |
---|
157 | wavelength = self.__getWavelength(frequency, float(wavelength)) |
---|
158 | |
---|
159 | return self._data[frequency][wavelength] |
---|
160 | else: |
---|
161 | return self.PDInfo(None) |
---|
162 | def getFocusPos(self): |
---|
163 | return self._focusPos |
---|
164 | |
---|
165 | def category(self): |
---|
166 | return "Diffraction;PythonAlgorithms" |
---|
167 | |
---|
168 | def name(self): |
---|
169 | return "SNSPowderReduction2" |
---|
170 | |
---|
171 | def PyInit(self): |
---|
172 | sns = mtd.getSettings().facility("SNS") |
---|
173 | instruments = [] |
---|
174 | for instr in sns.instruments(): |
---|
175 | for tech in instr.techniques(): |
---|
176 | if "Neutron Diffraction" == str(tech): |
---|
177 | instruments.append(instr.shortName()) |
---|
178 | break |
---|
179 | self.declareProperty("Instrument", "PG3", |
---|
180 | Validator=ListValidator(instruments)) |
---|
181 | #types = ["Event preNeXus", "Event NeXus"] |
---|
182 | #self.declareProperty("FileType", "Event NeXus", |
---|
183 | # Validator=ListValidator(types)) |
---|
184 | self.declareListProperty("RunNumber", [0], Validator=ArrayBoundedValidator(Lower=0)) |
---|
185 | extensions = [ "_histo.nxs", "_event.nxs", "_runinfo.xml"] |
---|
186 | self.declareProperty("Extension", "_event.nxs", |
---|
187 | Validator=ListValidator(extensions)) |
---|
188 | self.declareProperty("PreserveEvents", True, |
---|
189 | Description="Argument to supply to algorithms that can change from events to histograms.") |
---|
190 | self.declareProperty("Sum", False, |
---|
191 | Description="Sum the runs. Does nothing for characterization runs") |
---|
192 | self.declareProperty("PushDataPositive", "None", |
---|
193 | Description="Add a constant to the data that makes it positive over the whole range.", |
---|
194 | Validator=ListValidator(["None", "ResetToZero", "AddMinimum"])) |
---|
195 | self.declareProperty("BackgroundNumber", 0, Validator=BoundedValidator(Lower=-1), |
---|
196 | Description="If specified overrides value in CharacterizationRunsFile If -1 turns off correction.") |
---|
197 | self.declareProperty("VanadiumNumber", 0, Validator=BoundedValidator(Lower=-1), |
---|
198 | Description="If specified overrides value in CharacterizationRunsFile. If -1 turns off correction.") |
---|
199 | self.declareProperty("VanadiumBackgroundNumber", 0, Validator=BoundedValidator(Lower=-1)) |
---|
200 | self.declareFileProperty("CalibrationFile", "", FileAction.Load, |
---|
201 | [".cal"]) |
---|
202 | self.declareFileProperty("CharacterizationRunsFile", "", FileAction.OptionalLoad, |
---|
203 | ['.txt'], |
---|
204 | Description="File with characterization runs denoted") |
---|
205 | self.declareProperty("UnwrapRef", 0., |
---|
206 | Description="Reference total flight path for frame unwrapping. Zero skips the correction") |
---|
207 | self.declareProperty("LowResRef", 0., |
---|
208 | Description="Reference DIFC for resolution removal. Zero skips the correction") |
---|
209 | self.declareProperty("CropWavelengthMin", 0., |
---|
210 | Description="Crop the data at this minimum wavelength. Overrides LowResRef.") |
---|
211 | self.declareProperty("RemovePromptPulseWidth", 0.0, |
---|
212 | Description="Width of events (in microseconds) near the prompt pulse to remove. 0 disables") |
---|
213 | self.declareProperty("FilterByTimeMin", 0., |
---|
214 | Description="Relative time to start filtering by in seconds. Applies only to sample.") |
---|
215 | self.declareProperty("FilterByTimeMax", 0., |
---|
216 | Description="Relative time to stop filtering by in seconds. Applies only to sample.") |
---|
217 | self.declareProperty("MaxChunkSize", 0.0, Description="Specify maximum Gbytes of file to read in one chunk. Default is whole file.") |
---|
218 | self.declareProperty("FilterCharacterizations", False, |
---|
219 | Description="Filter the characterization runs using above parameters. This only works for event files.") |
---|
220 | self.declareListProperty("Binning", [0.,0.,0.], |
---|
221 | Description="Positive is linear bins, negative is logorithmic") |
---|
222 | self.declareProperty("BinInDspace", True, |
---|
223 | Description="If all three bin parameters a specified, whether they are in dspace (true) or time-of-flight (false)") |
---|
224 | self.declareProperty("StripVanadiumPeaks", True, |
---|
225 | Description="Subtract fitted vanadium peaks from the known positions.") |
---|
226 | self.declareProperty("VanadiumFWHM", 7, Description="Default=7") |
---|
227 | self.declareProperty("VanadiumPeakTol", 0.05, |
---|
228 | Description="How far from the ideal position a vanadium peak can be during StripVanadiumPeaks. Default=0.05, negative turns off") |
---|
229 | self.declareProperty("VanadiumSmoothParams", "20,2", Description="Default=20,2") |
---|
230 | self.declareProperty("FilterBadPulses", True, Description="Filter out events measured while proton charge is more than 5% below average") |
---|
231 | outfiletypes = ['gsas', 'fullprof', 'gsas and fullprof', 'gsas and fullprof and pdfgetn'] |
---|
232 | self.declareProperty("FilterByLogValue", "", Description="Name of log value to filter by") |
---|
233 | self.declareProperty("FilterMinimumValue", 0.0, Description="Minimum log value for which to keep events.") |
---|
234 | self.declareProperty("FilterMaximumValue", 0.0, Description="Maximum log value for which to keep events.") |
---|
235 | self.declareProperty("SaveAs", "gsas", ListValidator(outfiletypes)) |
---|
236 | self.declareFileProperty("OutputDirectory", "", FileAction.Directory) |
---|
237 | self.declareProperty("NormalizeByCurrent", True, Description="Normalized by Current") |
---|
238 | self.declareProperty("FinalDataUnits", "dSpacing", ListValidator(["dSpacing","MomentumTransfer"])) |
---|
239 | |
---|
240 | def _loadData(self, runnumber, extension, filterWall=None, **chunk): |
---|
241 | if runnumber is None or runnumber <= 0: |
---|
242 | return None |
---|
243 | |
---|
244 | name = "%s_%d" % (self._instrument, runnumber) |
---|
245 | filename = name + extension |
---|
246 | # EMPTY_INT() from C++ |
---|
247 | if int(chunk["ChunkNumber"]) < 2147483647: |
---|
248 | name += "_%02d" % (int(chunk["ChunkNumber"])) |
---|
249 | else: |
---|
250 | name += "_%02d" % 0 |
---|
251 | |
---|
252 | if extension.endswith("_event.nxs"): |
---|
253 | chunk["Precount"] = True |
---|
254 | if filterWall is not None: |
---|
255 | if filterWall[0] > 0.: |
---|
256 | chunk["FilterByTimeStart"] = filterWall[0] |
---|
257 | if filterWall[1] > 0.: |
---|
258 | chunk["FilterByTimeStop"] = filterWall[1] |
---|
259 | elif extension.endswith("_histo.nxs"): |
---|
260 | chunk = {} |
---|
261 | |
---|
262 | alg = Load(Filename=filename, OutputWorkspace=name, **chunk) |
---|
263 | return alg.workspace() |
---|
264 | |
---|
265 | def _getStrategy(self, runnumber, extension): |
---|
266 | # generate the workspace name |
---|
267 | wksp = "%s_%d" % (self._instrument, runnumber) |
---|
268 | strategy = [] |
---|
269 | Chunks = DetermineChunking(Filename=wksp+extension,MaxChunkSize=self._chunks,OutputWorkspace='Chunks').workspace() |
---|
270 | for row in Chunks: strategy.append(row) |
---|
271 | |
---|
272 | return strategy |
---|
273 | |
---|
274 | def _focusChunks(self, runnumber, extension, filterWall, calib, filterLogs=None, preserveEvents=True, |
---|
275 | normByCurrent=True, filterBadPulsesOverride=True): |
---|
276 | # generate the workspace name |
---|
277 | wksp = "%s_%d" % (self._instrument, runnumber) |
---|
278 | strategy = self._getStrategy(runnumber, extension) |
---|
279 | firstChunk = True |
---|
280 | for chunk in strategy: |
---|
281 | if "ChunkNumber" in chunk: |
---|
282 | self.log().information("Working on chunk %d of %d" % (chunk["ChunkNumber"], chunk["TotalChunks"])) |
---|
283 | temp = self._loadData(runnumber, extension, filterWall, **chunk) |
---|
284 | if self._info is None: |
---|
285 | self._info = self._getinfo(temp) |
---|
286 | temp = self._focus(temp, calib, self._info, filterLogs, preserveEvents, normByCurrent, filterBadPulsesOverride) |
---|
287 | if firstChunk: |
---|
288 | alg = RenameWorkspace(InputWorkspace=temp, OutputWorkspace=wksp) |
---|
289 | wksp = alg['OutputWorkspace'] |
---|
290 | firstChunk = False |
---|
291 | else: |
---|
292 | wksp += temp |
---|
293 | DeleteWorkspace(temp) |
---|
294 | # Sum workspaces for all mpi tasks |
---|
295 | if HAVE_MPI: |
---|
296 | alg = GatherWorkspaces(InputWorkspace=wksp, PreserveEvents=preserveEvents, AccumulationMethod="Add", OutputWorkspace=wksp) |
---|
297 | wksp = alg['OutputWorkspace'] |
---|
298 | if self._chunks > 0: |
---|
299 | # When chunks are added, proton charge is summed for all chunks |
---|
300 | wksp.getRun().integrateProtonCharge() |
---|
301 | mtd.deleteWorkspace('Chunks') |
---|
302 | if preserveEvents and not "histo" in self.getProperty("Extension"): |
---|
303 | CompressEvents(InputWorkspace=wksp, OutputWorkspace=wksp, Tolerance=COMPRESS_TOL_TOF) # 100ns |
---|
304 | if normByCurrent: |
---|
305 | NormaliseByCurrent(InputWorkspace=wksp, OutputWorkspace=wksp) |
---|
306 | wksp.getRun()['gsas_monitor'] = 1 |
---|
307 | |
---|
308 | self._save(wksp, self._info, False, True) |
---|
309 | self.log().information("Done focussing data") |
---|
310 | |
---|
311 | return wksp |
---|
312 | |
---|
313 | def _focus(self, wksp, calib, info, filterLogs=None, preserveEvents=True, |
---|
314 | normByCurrent=True, filterBadPulsesOverride=True): |
---|
315 | if wksp is None: |
---|
316 | return None |
---|
317 | |
---|
318 | # load the calibration file if the workspaces aren't already in memory |
---|
319 | if (not mtd.workspaceExists(self._instrument + "_offsets")) \ |
---|
320 | or (not mtd.workspaceExists(self._instrument + "_mask")) \ |
---|
321 | or (not mtd.workspaceExists(self._instrument + "_group")): |
---|
322 | whichones = {} |
---|
323 | whichones['MakeGroupingWorkspace'] = (not mtd.workspaceExists(self._instrument + "_group")) |
---|
324 | whichones['MakeOffsetsWorkspace'] = (not mtd.workspaceExists(self._instrument + "_offsets")) |
---|
325 | whichones['MakeMaskWorkspace'] = (not mtd.workspaceExists(self._instrument + "_mask")) |
---|
326 | LoadCalFile(InputWorkspace=wksp, CalFileName=calib, WorkspaceName=self._instrument, |
---|
327 | **whichones) |
---|
328 | |
---|
329 | if not "histo" in self.getProperty("Extension"): |
---|
330 | # take care of filtering events |
---|
331 | if self._filterBadPulses and filterBadPulsesOverride: |
---|
332 | FilterBadPulses(InputWorkspace=wksp, OutputWorkspace=wksp) |
---|
333 | if self._removePromptPulseWidth > 0.: |
---|
334 | RemovePromptPulse(InputWorkspace=wksp, OutputWorkspace=wksp, Width= self._removePromptPulseWidth) |
---|
335 | if filterLogs is not None: |
---|
336 | try: |
---|
337 | logparam = wksp.getRun()[filterLogs[0]] |
---|
338 | if logparam is not None: |
---|
339 | FilterByLogValue(InputWorkspace=wksp, OutputWorkspace=wksp, LogName=filterLogs[0], |
---|
340 | MinimumValue=filterLogs[1], MaximumValue=filterLogs[2]) |
---|
341 | except KeyError, e: |
---|
342 | raise RuntimeError("Failed to find log '%s' in workspace '%s'" \ |
---|
343 | % (filterLogs[0], str(wksp))) |
---|
344 | CompressEvents(InputWorkspace=wksp, OutputWorkspace=wksp, Tolerance=COMPRESS_TOL_TOF) # 100ns |
---|
345 | SortEvents(wksp) |
---|
346 | # determine some bits about d-space and binning |
---|
347 | if len(self._binning) == 3: |
---|
348 | info.has_dspace = self._bin_in_dspace |
---|
349 | if info.has_dspace: |
---|
350 | if len(self._binning) == 3: |
---|
351 | binning = self._binning |
---|
352 | else: |
---|
353 | binning = [info.dmin, self._binning[0], info.dmax] |
---|
354 | self.log().information("d-Spacing Binning: " + str(binning)) |
---|
355 | else: |
---|
356 | if len(self._binning) == 3: |
---|
357 | binning = self._binning |
---|
358 | else: |
---|
359 | binning = [info.tmin, self._binning[0], info.tmax] |
---|
360 | |
---|
361 | # trim off the ends |
---|
362 | cropkwargs = {} |
---|
363 | if info.tmin > 0.: |
---|
364 | cropkwargs["XMin"] = info.tmin |
---|
365 | if info.tmax > 0.: |
---|
366 | cropkwargs["XMax"] = info.tmax |
---|
367 | if not info.has_dspace: |
---|
368 | cropkwargs["XMin"] = binning[0] |
---|
369 | cropkwargs["XMax"] = binning[-1] |
---|
370 | if len(cropkwargs) > 0: |
---|
371 | CropWorkspace(InputWorkspace=wksp, OutputWorkspace=wksp, **cropkwargs) |
---|
372 | MaskDetectors(Workspace=wksp, MaskedWorkspace=self._instrument + "_mask") |
---|
373 | if info.has_dspace: |
---|
374 | if binning[0] > 100: |
---|
375 | self.log().warning("Binning data oddly for d-spacing: %s" % str(binning)) |
---|
376 | else: |
---|
377 | Rebin(InputWorkspace=wksp, OutputWorkspace=wksp, Params=binning) |
---|
378 | AlignDetectors(InputWorkspace=wksp, OutputWorkspace=wksp, OffsetsWorkspace=self._instrument + "_offsets") |
---|
379 | LRef = self.getProperty("UnwrapRef") |
---|
380 | DIFCref = self.getProperty("LowResRef") |
---|
381 | wavelengthMin = self.getProperty("CropWavelengthMin") |
---|
382 | if (LRef > 0.) or (DIFCref > 0.) or (wavelengthMin>0): # super special Jason stuff |
---|
383 | kwargs = {} |
---|
384 | try: |
---|
385 | if info.tmin > 0: |
---|
386 | kwargs["Tmin"] = info.tmin |
---|
387 | if info.tmax > info.tmin: |
---|
388 | kwargs["Tmax"] = info.tmax |
---|
389 | except: |
---|
390 | pass |
---|
391 | ConvertUnits(InputWorkspace=wksp, OutputWorkspace=wksp, Target="TOF") # corrections only work in TOF for now |
---|
392 | if LRef > 0.: |
---|
393 | UnwrapSNS(InputWorkspace=wksp, OutputWorkspace=wksp, LRef=LRef, **kwargs) |
---|
394 | if DIFCref > 0. or wavelengthMin > 0.: |
---|
395 | if kwargs.has_key("Tmax"): |
---|
396 | del kwargs["Tmax"] |
---|
397 | if wavelengthMin > 0.: |
---|
398 | RemoveLowResTOF(InputWorkspace=wksp, OutputWorkspace=wksp, MinWavelength=wavelengthMin, **kwargs) |
---|
399 | else: |
---|
400 | RemoveLowResTOF(InputWorkspace=wksp, OutputWorkspace=wksp, ReferenceDIFC=DIFCref, |
---|
401 | K=3.22, **kwargs) |
---|
402 | ConvertUnits(InputWorkspace=wksp, OutputWorkspace=wksp, Target="dSpacing") # put it back to the original units |
---|
403 | if info.has_dspace: |
---|
404 | self.log().information("d-Spacing binning: " + str(binning)) |
---|
405 | Rebin(InputWorkspace=wksp, OutputWorkspace=wksp, Params=binning) |
---|
406 | if not "histo" in self.getProperty("Extension"): |
---|
407 | SortEvents(InputWorkspace=wksp) |
---|
408 | DiffractionFocussing(InputWorkspace=wksp, OutputWorkspace=wksp, GroupingWorkspace=self._instrument + "_group", |
---|
409 | PreserveEvents=preserveEvents) |
---|
410 | if not "histo" in self.getProperty("Extension") and preserveEvents: |
---|
411 | SortEvents(InputWorkspace=wksp) |
---|
412 | |
---|
413 | focusPos = self._config.getFocusPos() |
---|
414 | self.log().notice("FOCUS:" + str(focusPos)) |
---|
415 | if not focusPos is None: |
---|
416 | EditInstrumentGeometry(Workspace=wksp, NewInstrument=False, **focusPos) |
---|
417 | if (self._config.iparmFile is not None) and (len(self._config.iparmFile) > 0): |
---|
418 | wksp.getRun()['iparm_file'] = self._config.iparmFile |
---|
419 | |
---|
420 | ConvertUnits(InputWorkspace=wksp, OutputWorkspace=wksp, Target="TOF") |
---|
421 | Rebin(InputWorkspace=wksp, OutputWorkspace=wksp, Params=binning[1]) # reset bin width |
---|
422 | return wksp |
---|
423 | |
---|
424 | def _getinfo(self, wksp): |
---|
425 | logs = wksp.getRun() |
---|
426 | # get the frequency |
---|
427 | frequency = None |
---|
428 | if "SpeedRequest1" in logs.keys(): |
---|
429 | frequency = logs['SpeedRequest1'] |
---|
430 | else: |
---|
431 | self.log().information("'SpeedRequest1' is not specified in logs") |
---|
432 | if "frequency" in logs.keys(): |
---|
433 | frequency = logs['frequency'] |
---|
434 | else: |
---|
435 | self.log().information("'frequency' is not specified in logs") |
---|
436 | return self._config.getInfo(None, None) |
---|
437 | if frequency.units != "Hz": |
---|
438 | raise RuntimeError("Only know how to deal with frequency in Hz, not %s" % frequency.units) |
---|
439 | frequency = frequency.getStatistics().mean |
---|
440 | |
---|
441 | if not "LambdaRequest" in logs.keys(): |
---|
442 | self.log().information("'LambdaRequest' is not in the datafile") |
---|
443 | return self._config.getInfo(None, None) |
---|
444 | wavelength = logs['LambdaRequest'] |
---|
445 | if wavelength.units != "Angstrom": |
---|
446 | raise RuntimeError("Only know how to deal with LambdaRequest in Angstrom, not $s" % wavelength) |
---|
447 | wavelength = wavelength.getStatistics().mean |
---|
448 | |
---|
449 | self.log().information("frequency: " + str(frequency) + "Hz center wavelength:" + str(wavelength) + "Angstrom") |
---|
450 | return self._config.getInfo(frequency, wavelength) |
---|
451 | |
---|
452 | def _save(self, wksp, info, normalized, pdfgetn): |
---|
453 | filename = os.path.join(self._outDir, str(wksp)) |
---|
454 | if pdfgetn: |
---|
455 | if "pdfgetn" in self._outTypes: |
---|
456 | pdfwksp = str(wksp)+"_norm" |
---|
457 | SetUncertainties(InputWorkspace=wksp, OutputWorkspace=pdfwksp, SetError="sqrt") |
---|
458 | SaveGSS(InputWorkspace=pdfwksp, Filename=filename+".getn", SplitFiles="False", Append=False, |
---|
459 | MultiplyByBinWidth=False, Bank=info.bank, Format="SLOG", ExtendedHeader=True) |
---|
460 | mtd.deleteWorkspace(pdfwksp) |
---|
461 | return # don't do the other bits of saving |
---|
462 | if "gsas" in self._outTypes: |
---|
463 | SaveGSS(InputWorkspace=wksp, Filename=filename+".gsa", SplitFiles="False", Append=False, |
---|
464 | MultiplyByBinWidth=normalized, Bank=info.bank, Format="SLOG", ExtendedHeader=True) |
---|
465 | if "fullprof" in self._outTypes: |
---|
466 | SaveFocusedXYE(InputWorkspace=wksp, StartAtBankNumber=info.bank, Filename=filename+".dat") |
---|
467 | |
---|
468 | # always save python script |
---|
469 | GeneratePythonScript(InputWorkspace=wksp, Filename=filename+".py") |
---|
470 | |
---|
471 | def PyExec(self): |
---|
472 | # get generic information |
---|
473 | SUFFIX = self.getProperty("Extension") |
---|
474 | self._config = self.PDConfigFile(self.getProperty("CharacterizationRunsFile")) |
---|
475 | self._binning = self.getProperty("Binning") |
---|
476 | if len(self._binning) != 1 and len(self._binning) != 3: |
---|
477 | raise RuntimeError("Can only specify (width) or (start,width,stop) for binning. Found %d values." % len(self._binning)) |
---|
478 | if len(self._binning) == 3: |
---|
479 | if self._binning[0] == 0. and self._binning[1] == 0. and self._binning[2] == 0.: |
---|
480 | raise RuntimeError("Failed to specify the binning") |
---|
481 | self._bin_in_dspace = self.getProperty("BinInDspace") |
---|
482 | self._instrument = self.getProperty("Instrument") |
---|
483 | mtd.settings['default.facility'] = "SNS" |
---|
484 | mtd.settings['default.instrument'] = self._instrument |
---|
485 | # self._timeMin = self.getProperty("FilterByTimeMin") |
---|
486 | # self._timeMax = self.getProperty("FilterByTimeMax") |
---|
487 | self._filterBadPulses = self.getProperty("FilterBadPulses") |
---|
488 | self._removePromptPulseWidth = self.getProperty("RemovePromptPulseWidth") |
---|
489 | filterLogs = self.getProperty("FilterByLogValue") |
---|
490 | if len(filterLogs.strip()) <= 0: |
---|
491 | filterLogs = None |
---|
492 | else: |
---|
493 | filterLogs = [filterLogs, |
---|
494 | self.getProperty("FilterMinimumValue"), self.getProperty("FilterMaximumValue")] |
---|
495 | self._vanPeakFWHM = self.getProperty("VanadiumFWHM") |
---|
496 | self._vanSmoothing = self.getProperty("VanadiumSmoothParams") |
---|
497 | calib = self.getProperty("CalibrationFile") |
---|
498 | self._outDir = self.getProperty("OutputDirectory") |
---|
499 | self._outTypes = self.getProperty("SaveAs") |
---|
500 | samRuns = self.getProperty("RunNumber") |
---|
501 | filterWall = (self.getProperty("FilterByTimeMin"), self.getProperty("FilterByTimeMax")) |
---|
502 | preserveEvents = self.getProperty("PreserveEvents") |
---|
503 | normbycurrent = self.getProperty("NormalizeByCurrent") |
---|
504 | self._info = None |
---|
505 | self._chunks = self.getProperty("MaxChunkSize") |
---|
506 | |
---|
507 | workspacelist = [] # all data workspaces that will be converted to d-spacing in the end |
---|
508 | |
---|
509 | if self.getProperty("Sum"): |
---|
510 | samRun = None |
---|
511 | info = None |
---|
512 | for temp in samRuns: |
---|
513 | temp = self._focusChunks(temp, SUFFIX, filterWall, calib, filterLogs, |
---|
514 | preserveEvents=preserveEvents, normByCurrent=normbycurrent) |
---|
515 | tempinfo = self._getinfo(temp) |
---|
516 | if samRun is None: |
---|
517 | samRun = temp |
---|
518 | info = tempinfo |
---|
519 | else: |
---|
520 | if (tempinfo.freq is not None) and (info.freq is not None) \ |
---|
521 | and (abs(tempinfo.freq - info.freq)/info.freq > .05): |
---|
522 | raise RuntimeError("Cannot add incompatible frequencies (%f!=%f)" \ |
---|
523 | % (tempinfo.freq, info.freq)) |
---|
524 | if (tempinfo.wl is not None) and (info.wl is not None) \ |
---|
525 | and abs(tempinfo.wl - info.wl)/info.freq > .05: |
---|
526 | raise RuntimeError("Cannot add incompatible wavelengths (%f != %f)" \ |
---|
527 | % (tempinfo.wl, info.wl)) |
---|
528 | Plus(samRun, temp, samRun) |
---|
529 | if not "histo" in self.getProperty("Extension"): |
---|
530 | CompressEvents(InputWorkspace=samRun, OutputWorkspace=samRun, |
---|
531 | Tolerance=COMPRESS_TOL_TOF) # 10ns |
---|
532 | mtd.deleteWorkspace(str(temp)) |
---|
533 | samRun /= float(len(samRuns)) |
---|
534 | samRuns = [samRun] |
---|
535 | workspacelist.append(str(samRun)) |
---|
536 | |
---|
537 | for samRun in samRuns: |
---|
538 | # first round of processing the sample |
---|
539 | if not self.getProperty("Sum"): |
---|
540 | self._info = None |
---|
541 | samRun = self._focusChunks(samRun, SUFFIX, filterWall, calib, filterLogs, |
---|
542 | preserveEvents=preserveEvents, normByCurrent=normbycurrent) |
---|
543 | workspacelist.append(str(samRun)) |
---|
544 | |
---|
545 | # process the container |
---|
546 | canRun = self.getProperty("BackgroundNumber") |
---|
547 | if canRun == 0: # use the version in the info |
---|
548 | canRun = self._info.can |
---|
549 | elif canRun < 0: # turn off the correction |
---|
550 | canRun = 0 |
---|
551 | if canRun > 0: |
---|
552 | if mtd.workspaceExists("%s_%d" % (self._instrument, canRun)): |
---|
553 | canRun = mtd["%s_%d" % (self._instrument, canRun)] |
---|
554 | else: |
---|
555 | if self.getProperty("FilterCharacterizations"): |
---|
556 | canRun = self._focusChunks(canRun, SUFFIX, filterWall, calib, |
---|
557 | preserveEvents=preserveEvents) |
---|
558 | else: |
---|
559 | canRun = self._focusChunks(canRun, SUFFIX, (0., 0.), calib, |
---|
560 | preserveEvents=preserveEvents) |
---|
561 | if HAVE_MPI: |
---|
562 | if mpi.world.rank == 0: |
---|
563 | ConvertUnits(InputWorkspace=canRun, OutputWorkspace=canRun, Target="TOF") |
---|
564 | workspacelist.append(str(canRun)) |
---|
565 | else: |
---|
566 | ConvertUnits(InputWorkspace=canRun, OutputWorkspace=canRun, Target="TOF") |
---|
567 | workspacelist.append(str(canRun)) |
---|
568 | else: |
---|
569 | canRun = None |
---|
570 | |
---|
571 | # process the vanadium run |
---|
572 | vanRun = self.getProperty("VanadiumNumber") |
---|
573 | if vanRun == 0: # use the version in the info |
---|
574 | vanRun = self._info.van |
---|
575 | elif vanRun < 0: # turn off the correction |
---|
576 | vanRun = 0 |
---|
577 | if vanRun > 0: |
---|
578 | if mtd.workspaceExists("%s_%d" % (self._instrument, vanRun)): |
---|
579 | vanRun = mtd["%s_%d" % (self._instrument, vanRun)] |
---|
580 | else: |
---|
581 | vnoiseRun = self._info.vnoise # noise run for the vanadium |
---|
582 | if self.getProperty("FilterCharacterizations"): |
---|
583 | vanRun = self._focusChunks(vanRun, SUFFIX, filterWall, calib, |
---|
584 | preserveEvents=False, normByCurrent = (vnoiseRun <= 0)) |
---|
585 | else: |
---|
586 | vanRun = self._focusChunks(vanRun, SUFFIX, (0., 0.), calib, |
---|
587 | preserveEvents=False, normByCurrent = (vnoiseRun <= 0)) |
---|
588 | |
---|
589 | if (vnoiseRun > 0): |
---|
590 | if self.getProperty("FilterCharacterizations"): |
---|
591 | vnoiseRun = self._focusChunks(vnoiseRun, SUFFIX, filterWall, calib, |
---|
592 | preserveEvents=False, normByCurrent = False, filterBadPulsesOverride=False) |
---|
593 | else: |
---|
594 | vnoiseRun = self._focusChunks(vnoiseRun, SUFFIX, (0., 0.), calib, |
---|
595 | preserveEvents=False, normByCurrent = False, filterBadPulsesOverride=False) |
---|
596 | if HAVE_MPI: |
---|
597 | if mpi.world.rank == 0: |
---|
598 | ConvertUnits(InputWorkspace=vnoiseRun, OutputWorkspace=vnoiseRun, Target="TOF") |
---|
599 | FFTSmooth(InputWorkspace=vnoiseRun, OutputWorkspace=vnoiseRun, Filter="Butterworth", |
---|
600 | Params=self._vanSmoothing,IgnoreXBins=True,AllSpectra=True) |
---|
601 | try: |
---|
602 | vanDuration = vanRun.getRun().get('duration') |
---|
603 | vanDuration = vanDuration.value |
---|
604 | except: |
---|
605 | vanDuration = 1. |
---|
606 | try: |
---|
607 | vbackDuration = vnoiseRun.getRun().get('duration') |
---|
608 | vbackDuration = vbackDuration.value |
---|
609 | except: |
---|
610 | vbackDuration = 1. |
---|
611 | vnoiseRun *= (vanDuration/vbackDuration) |
---|
612 | vanRun -= vnoiseRun |
---|
613 | NormaliseByCurrent(InputWorkspace=vanRun, OutputWorkspace=vanRun) |
---|
614 | workspacelist.append(str(vnoiseRun)) |
---|
615 | else: |
---|
616 | ConvertUnits(InputWorkspace=vnoiseRun, OutputWorkspace=vnoiseRun, Target="TOF") |
---|
617 | FFTSmooth(InputWorkspace=vnoiseRun, OutputWorkspace=vnoiseRun, Filter="Butterworth", |
---|
618 | Params=self._vanSmoothing,IgnoreXBins=True,AllSpectra=True) |
---|
619 | try: |
---|
620 | vanDuration = vanRun.getRun().get('duration') |
---|
621 | vanDuration = vanDuration.value |
---|
622 | except: |
---|
623 | vanDuration = 1. |
---|
624 | try: |
---|
625 | vbackDuration = vnoiseRun.getRun().get('duration') |
---|
626 | vbackDuration = vbackDuration.value |
---|
627 | except: |
---|
628 | vbackDuration = 1. |
---|
629 | vnoiseRun *= (vanDuration/vbackDuration) |
---|
630 | vanRun -= vnoiseRun |
---|
631 | NormaliseByCurrent(InputWorkspace=vanRun, OutputWorkspace=vanRun) |
---|
632 | workspacelist.append(str(vnoiseRun)) |
---|
633 | else: |
---|
634 | vnoiseRun = None |
---|
635 | |
---|
636 | vbackRun = self.getProperty("VanadiumBackgroundNumber") |
---|
637 | if vbackRun > 0: |
---|
638 | if mtd.workspaceExists("%s_%d" % (self._instrument, vbackRun)): |
---|
639 | vbackRun = mtd["%s_%d" % (self._instrument, vbackRun)] |
---|
640 | else: |
---|
641 | if self.getProperty("FilterCharacterizations"): |
---|
642 | vbackRun = self._focusChunks(vbackRun, SUFFIX, filterWall, calib, |
---|
643 | preserveEvents=False) |
---|
644 | else: |
---|
645 | vbackRun = self._focusChunks(vbackRun, SUFFIX, (0., 0.), calib, |
---|
646 | preserveEvents=False) |
---|
647 | vanRun -= vbackRun |
---|
648 | workspacelist.append(str(vbackRun)) |
---|
649 | |
---|
650 | if HAVE_MPI: |
---|
651 | if mpi.world.rank > 0: |
---|
652 | return |
---|
653 | if self.getProperty("StripVanadiumPeaks"): |
---|
654 | ConvertUnits(InputWorkspace=vanRun, OutputWorkspace=vanRun, Target="dSpacing") |
---|
655 | StripVanadiumPeaks(InputWorkspace=vanRun, OutputWorkspace=vanRun, FWHM=self._vanPeakFWHM, |
---|
656 | PeakPositionTolerance=self.getProperty("VanadiumPeakTol"), |
---|
657 | BackgroundType="Quadratic", HighBackground=True) |
---|
658 | ConvertUnits(InputWorkspace=vanRun, OutputWorkspace=vanRun, Target="TOF") |
---|
659 | FFTSmooth(InputWorkspace=vanRun, OutputWorkspace=vanRun, Filter="Butterworth", |
---|
660 | Params=self._vanSmoothing,IgnoreXBins=True,AllSpectra=True) |
---|
661 | MultipleScatteringCylinderAbsorption(InputWorkspace=vanRun, OutputWorkspace=vanRun, # numbers for vanadium |
---|
662 | AttenuationXSection=2.8, ScatteringXSection=5.1, |
---|
663 | SampleNumberDensity=0.0721, CylinderSampleRadius=.3175) |
---|
664 | SetUncertainties(InputWorkspace=vanRun, OutputWorkspace=vanRun) |
---|
665 | if HAVE_MPI: |
---|
666 | if mpi.world.rank > 0: |
---|
667 | return |
---|
668 | ConvertUnits(InputWorkspace=vanRun, OutputWorkspace=vanRun, Target="TOF") |
---|
669 | workspacelist.append(str(vanRun)) |
---|
670 | else: |
---|
671 | vanRun = None |
---|
672 | |
---|
673 | # the final bit of math |
---|
674 | if canRun is not None: |
---|
675 | samRun -= canRun |
---|
676 | if not "histo" in self.getProperty("Extension") and preserveEvents: |
---|
677 | CompressEvents(InputWorkspace=samRun, OutputWorkspace=samRun, |
---|
678 | Tolerance=COMPRESS_TOL_TOF) # 10ns |
---|
679 | canRun = str(canRun) |
---|
680 | if vanRun is not None: |
---|
681 | samRun /= vanRun |
---|
682 | normalized = True |
---|
683 | samRun.getRun()['van_number'] = vanRun.getRun()['run_number'].value |
---|
684 | vanRun = str(vanRun) |
---|
685 | else: |
---|
686 | normalized = False |
---|
687 | |
---|
688 | if not "histo" in self.getProperty("Extension") and preserveEvents and HAVE_MPI is False: |
---|
689 | CompressEvents(InputWorkspace=samRun, OutputWorkspace=samRun, |
---|
690 | Tolerance=COMPRESS_TOL_TOF) # 5ns |
---|
691 | |
---|
692 | # make sure there are no negative values - gsas hates them |
---|
693 | if self.getProperty("PushDataPositive") != "None": |
---|
694 | addMin = (self.getProperty("PushDataPositive") == "AddMinimum") |
---|
695 | ResetNegatives(InputWorkspace=samRun, OutputWorkspace=samRun, AddMinimum=addMin, ResetValue=0.) |
---|
696 | |
---|
697 | # write out the files |
---|
698 | if HAVE_MPI: |
---|
699 | if mpi.world.rank == 0: |
---|
700 | self._save(samRun, self._info, normalized, False) |
---|
701 | samRun = str(samRun) |
---|
702 | else: |
---|
703 | self._save(samRun, self._info, normalized, False) |
---|
704 | samRun = str(samRun) |
---|
705 | mtd.releaseFreeMemory() |
---|
706 | |
---|
707 | # convert everything into d-spacing |
---|
708 | workspacelist = set(workspacelist) # only do each workspace once |
---|
709 | if HAVE_MPI is False: |
---|
710 | for wksp in workspacelist: |
---|
711 | ConvertUnits(InputWorkspace=wksp, OutputWorkspace=wksp, Target=self.getProperty("FinalDataUnits")) |
---|
712 | |
---|
713 | mtd.registerPyAlgorithm(SNSPowderReduction2()) |
---|