Ticket #5743: SNSPowderReduction2.py

File SNSPowderReduction2.py, 36.1 KB (added by Vickie Lynch, 8 years ago)

SNSPowderReduction before changes

Line 
1"""*WIKI*
2
3
4*WIKI*"""
5
6from MantidFramework import *
7from mantidsimple import *
8from mantid.api import AlgorithmFactory
9import os
10
11all_algs = AlgorithmFactory.getRegisteredAlgorithms(True)
12if 'GatherWorkspaces' in all_algs:
13    HAVE_MPI = True
14    import boostmpi as mpi
15else:
16    HAVE_MPI = False
17
18COMPRESS_TOL_TOF = .01
19
20class 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
713mtd.registerPyAlgorithm(SNSPowderReduction2())