diff --git a/Code/Mantid/instrument/INTER_Definition.xml b/Code/Mantid/instrument/INTER_Definition.xml index 67f9e3842edb..ff9ce84c7f3d 100644 --- a/Code/Mantid/instrument/INTER_Definition.xml +++ b/Code/Mantid/instrument/INTER_Definition.xml @@ -41,7 +41,7 @@ - + diff --git a/Code/Mantid/scripts/CMakeLists.txt b/Code/Mantid/scripts/CMakeLists.txt index fea04ffa06dd..815236738bc9 100644 --- a/Code/Mantid/scripts/CMakeLists.txt +++ b/Code/Mantid/scripts/CMakeLists.txt @@ -1,11 +1,11 @@ set ( TEST_PY_FILES + test/ConvertToWavelengthTest.py test/ReducerTest.py test/SettingsTest.py test/DgreduceTest.py test/DirectEnergyConversionTest.py test/ReflectometryQuickAuxiliaryTest.py - test/ReflectometryQuickToLamTest.py test/ExtendedUnitCellTest.py test/SansIsisGuiSettings.py ) diff --git a/Code/Mantid/scripts/Reflectometry/isis_reflectometry/convert_to_wavelength.py b/Code/Mantid/scripts/Reflectometry/isis_reflectometry/convert_to_wavelength.py new file mode 100644 index 000000000000..7f5d4c9bb79f --- /dev/null +++ b/Code/Mantid/scripts/Reflectometry/isis_reflectometry/convert_to_wavelength.py @@ -0,0 +1,133 @@ +import mantid.simpleapi as msi +import mantid.api +from mantid.kernel import logger + +class ConvertToWavelength(object): + + # List of workspaces to process. + __ws_list = [] + + @classmethod + def sum_workspaces(cls, workspaces): + """ + Sum together all workspaces. return the result. + Returns: + Result of sum ( a workspace) + """ + return sum(workspaces) + + @classmethod + def to_workspace(cls, candidate): + workspace = None + if isinstance(candidate, mantid.api.MatrixWorkspace): + workspace = candidate + elif isinstance(candidate, str): + if mantid.api.AnalysisDataService.doesExist(candidate): + workspace = mantid.api.AnalysisDataService.retrieve(candidate) + else: + workspace = msi.Load(Filename=candidate) + else: + raise ValueError("Unknown source item %s" % candidate) + return workspace + + def __to_workspace_list(self, source_list): + temp=[] + for item in source_list: + temp.append(ConvertToWavelength.to_workspace(item)) + self.__ws_list = temp + + + def __init__(self, source): + """ + Constructor + Arguments: + list -- source workspace or workspaces. + + Convert inputs into a list of workspace objects. + """ + source_list = None + if not isinstance(source, list): + source_list = [source] + else: + source_list = source + self.__to_workspace_list(source_list) + + @classmethod + def crop_range(cls, ws, rng): + """ + Given a range of workspace indexes to keep (rng), Crop the workspace such that these are kept. + Arguments: + + ws : Workspace to crop spectrum from + rng : Tuple of range tuples. syntax may be (start, stop) or ((start0, stop0), (start1, stop1), ...) + + returns a new copy of the workspace with spectra cropped out. + """ + + in_rng = msi.CloneWorkspace(ws) + if not isinstance(rng, tuple): + raise ValueError("Elements must be tuples.") + + def is_actual_range(rng): + return len(rng) == 2 and isinstance(rng[0], int) and isinstance(rng[1], int) + + if is_actual_range(rng): + start, stop = rng[0], rng[1] + in_rng = msi.CropWorkspace(InputWorkspace=in_rng, StartWorkspaceIndex=start,EndWorkspaceIndex=stop) + else: + for subrng in rng: + in_rng = ConvertToWavelength.crop_range(ws, subrng) + + return in_rng + + def convert(self, wavelength_min, wavelength_max, detector_workspace_indexes, monitor_workspace_index, correct_monitor=False, bg_min=None, bg_max=None): + """ + Run the conversion + + Arguments: + + workspace_ids: Start and end ranges. Ids to be considered as workspaces. Nested list syntax supported + wavelength_min: min wavelength in x for monitor workspace + wavelength_max: max wavelength in x for detector workspace + detector_workspace_indexes: Tuple of workspace indexes (or tuple of tuple min, max ranges to keep) + monitor_workspace_index: The index of the monitor workspace + correct_monitor: Flag indicating that monitors should have a flat background correction applied + bg_min: x min background in wavelength + bg_max: x max background in wavelength + + Returns: + monitor_ws: A workspace of monitors + """ + # Sanity check inputs. + if(wavelength_min >= wavelength_max): + raise ValueError("Wavelength_min must be < wavelength_max min: %s, max: %s" % (wavelength_min, wavelength_max)) + + if correct_monitor and not all((bg_min, bg_max)): + raise ValueError("Either provide ALL, monitors_to_correct, bg_min, bg_max or none of them") + + if all((bg_min, bg_max)) and bg_min >= bg_max: + raise ValueError("Background min must be < Background max") + + sum = ConvertToWavelength.sum_workspaces(self.__ws_list) + sum_wavelength= msi.ConvertUnits(InputWorkspace=sum, Target="Wavelength", AlignBins='1') + + logger.debug("Monitor detector index %s" % str(monitor_workspace_index)) + + # Crop out the monitor workspace + monitor_ws = msi.CropWorkspace(InputWorkspace=sum_wavelength, StartWorkspaceIndex=monitor_workspace_index,EndWorkspaceIndex=monitor_workspace_index) + # Crop out the detector workspace then chop out the x-ranges of interest. + detector_ws = ConvertToWavelength.crop_range(sum_wavelength, detector_workspace_indexes) + + detector_ws = msi.CropWorkspace(InputWorkspace=detector_ws, XMin=wavelength_min, XMax=wavelength_max) + + # Apply a flat background + if correct_monitor and all((bg_min, bg_max)): + monitor_ws = msi.CalculateFlatBackground(InputWorkspace=monitor_ws,WorkspaceIndexList=0,StartX=bg_min, EndX=bg_max) + + msi.DeleteWorkspace(Workspace=sum_wavelength.getName()) + return (monitor_ws, detector_ws) + + + + + diff --git a/Code/Mantid/scripts/Reflectometry/isis_reflectometry/quick.py b/Code/Mantid/scripts/Reflectometry/isis_reflectometry/quick.py index 726061602984..8e9f9c125724 100644 --- a/Code/Mantid/scripts/Reflectometry/isis_reflectometry/quick.py +++ b/Code/Mantid/scripts/Reflectometry/isis_reflectometry/quick.py @@ -1,9 +1,9 @@ -''' SVN Info: The variables below will only get subsituted at svn checkout if - the repository is configured for variable subsitution. +''' SVN Info: The variables below will only get subsituted at svn checkout if + the repository is configured for variable subsitution. - $Id$ - $HeadURL$ -|=============================================================================|=======| + $Id$ + $HeadURL$ +|=============================================================================|=======| 1 80 ''' #these need to be moved into one NR folder or so @@ -13,468 +13,385 @@ #from mantidsimple import * # Old API from mantid.simpleapi import * # New API from mantid.api import WorkspaceGroup +from convert_to_wavelength import ConvertToWavelength import math import re def quick(run, theta=0, pointdet=1,roi=[0,0], db=[0,0], trans='', polcorr=0, usemon=-1,outputType='pd', debug=0): - ''' - call signature(s):: - - x=quick(RunNumber) - x=quick(RunNumber, roi=[0,0], db=[0,0], trans=0, outputType='pd') - x=quick(RunNumber,[1,10]) - x=quick(RunNumber,[1,10],[20,40]) - x=quick(RunNumber, trans=2568) - x=quick(RunNumber, trans='SomeSavedWorkspaceName') - - Reduces a ISIS raw or nexus file created on one of the reflectometers applying - only a minimum ammount of corrections. The data is left in terms of lambda. - - Required arguments - ========= ===================================================================== - RunNumber Either an ISIS run number when the paths are set up correctly or - the full path and filename if an ISIS raw file. - ========= ===================================================================== - - Optional keyword arguments: - ========= ===================================================================== - Keyword Description - ========= ===================================================================== - roi Region of interest marking the extent of the reflected beam. - default [0,0] - db Region of interest marking the extent of the direct beam. - default [0,0] - trans transmission run number or saved workspace. The default is 0 (No - transmission run). trans=-1 will supress the division of the - detector by the monitor. - polcorr polarisation correction, 0=no correction (unpolarised run) - usemon monitor to be used for normalisation (-1 is default from IDF) - outputType 'pd' = point detector (Default), 'md'= Multidetector Will use - this to build the equivalent of gd in the old matlab code but - keep all of the simple detector processing in one well organized - function. This should not be used by the average user. - ========= ===================================================================== - - Outputs: - ========= ===================================================================== - x Either a single mantid workspace or worspace group or an array - of them. - ========= ===================================================================== - - Working examples: - >>> # reduce a data set with the default parameters - >>> x=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003014.raw") - - >>> # reduce a data set with a transmission run - >>> t=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003010.raw") - >>> x=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003014.raw", trans=t) - - >>> # reduce a data set using the multidetector and output a single reflectivity - >>> # where the reflected beam is between channel 121 and 130. - >>> x=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003014.raw", [121,130]) + ''' + call signature(s):: + + x=quick(RunNumber) + x=quick(RunNumber, roi=[0,0], db=[0,0], trans=0, outputType='pd') + x=quick(RunNumber,[1,10]) + x=quick(RunNumber,[1,10],[20,40]) + x=quick(RunNumber, trans=2568) + x=quick(RunNumber, trans='SomeSavedWorkspaceName') + + Reduces a ISIS raw or nexus file created on one of the reflectometers applying + only a minimum ammount of corrections. The data is left in terms of lambda. + + Required arguments + ========= ===================================================================== + RunNumber Either an ISIS run number when the paths are set up correctly or + the full path and filename if an ISIS raw file. + ========= ===================================================================== + + Optional keyword arguments: + ========= ===================================================================== + Keyword Description + ========= ===================================================================== + roi Region of interest marking the extent of the reflected beam. + default [0,0] + db Region of interest marking the extent of the direct beam. + default [0,0] + trans transmission run number or saved workspace. The default is 0 (No + transmission run). trans=-1 will supress the division of the + detector by the monitor. + polcorr polarisation correction, 0=no correction (unpolarised run) + usemon monitor to be used for normalisation (-1 is default from IDF) + outputType 'pd' = point detector (Default), 'md'= Multidetector Will use + this to build the equivalent of gd in the old matlab code but + keep all of the simple detector processing in one well organized + function. This should not be used by the average user. + ========= ===================================================================== + + Outputs: + ========= ===================================================================== + x Either a single mantid workspace or worspace group or an array + of them. + ========= ===================================================================== + + Working examples: + >>> # reduce a data set with the default parameters + >>> x=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003014.raw") + + >>> # reduce a data set with a transmission run + >>> t=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003010.raw") + >>> x=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003014.raw", trans=t) + + >>> # reduce a data set using the multidetector and output a single reflectivity + >>> # where the reflected beam is between channel 121 and 130. + >>> x=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003014.raw", [121,130]) - Also see: pol + Also see: pol - ToDo: - 1) code for the transmisson DONE! - 2) Similar to the genie on polref add extraction from the multidetector - 3) need to make the variables stored in the frame work contain the run number. DONE! - - ''' + ToDo: + 1) code for the transmisson DONE! + 2) Similar to the genie on polref add extraction from the multidetector + 3) need to make the variables stored in the frame work contain the run number. DONE! + + ''' - ''' Notes for developers: + ''' Notes for developers: - Naming conventions for workspaces which live in the mantid framework are as follows: + Naming conventions for workspaces which live in the mantid framework are as follows: - It's nearly random. this needs to be fixed so that name clashes do not occur. - May try adding a pair of underscores to the front of the name. + It's nearly random. this needs to be fixed so that name clashes do not occur. + May try adding a pair of underscores to the front of the name. - ''' - - - [I0MonitorIndex, MultiDetectorStart, nHist] = toLam(run,'',pointdet=1) #creates wTof = "_W" + name - inst = groupGet("_W",'inst') - # Some beamline constants from IDF - intmin = inst.getNumberParameter('MonitorIntegralMin')[0] - intmax = inst.getNumberParameter('MonitorIntegralMax')[0] - print I0MonitorIndex - print nHist - if (nHist > 5 and not(pointdet)): - # Proccess Multi-Detector; assume MD goes to the end: - # if roi or db are given in the function then sum over the apropriate channels - print "This is a multidetector run." - try: - CropWorkspace(InputWorkspace="_D",OutputWorkspace="_DM",StartWorkspaceIndex=MultiDetectorStart) - RebinToWorkspace(WorkspaceToRebin="_M",WorkspaceToMatch="_DM",OutputWorkspace="_M_M") - CropWorkspace(InputWorkspace="_M_M",OutputWorkspace="_I0M",StartWorkspaceIndex=I0MonitorIndex) - RebinToWorkspace(WorkspaceToRebin="_I0M",WorkspaceToMatch="_DM",OutputWorkspace="_I0M") - Divide(LHSWorkspace="_DM",RHSWorkspace="_I0M",OutputWorkspace="IvsLam") - if (roi != [0,0]) : - SumSpectra(InputWorkspace="IvsLam",OutputWorkspace="DMR",StartWorkspaceIndex=roi[0], EndWorkspaceIndex=roi[1]) - ReflectedBeam=mtd['DMR'] - if (db != [0,0]) : - SumSpectra(InputWorkspace="_DM",OutputWorkspace="_DMD",StartWorkspaceIndex=db[0], EndWorkspaceIndex=db[1]) - DirectBeam=mtd['_DMD'] - except SystemExit: - print "Point-Detector only run." - RunNumber = groupGet('IvsLam','samp','run_number') - if (theta): - IvsQ = l2q(mtd['DMR'], 'linear-detector', theta) - else: - ConvertUnits(InputWorkspace='DMR',OutputWorkspace="IvsQ",Target="MomentumTransfer") - - # Single Detector processing------------------------------------------------------------- - else: - print "This is a Point-Detector run." - # handle transmission runs - # process the point detector reflectivity - RebinToWorkspace(WorkspaceToRebin="_M",WorkspaceToMatch="_DP",OutputWorkspace="_M_P") - CropWorkspace(InputWorkspace="_M_P",OutputWorkspace="_I0P",StartWorkspaceIndex=I0MonitorIndex,EndWorkspaceIndex=I0MonitorIndex) - Scale(InputWorkspace="_DP",OutputWorkspace="IvsLam",Factor=1) - #Divide(LHSWorkspace="_DP",RHSWorkspace="_I0P",OutputWorkspace="IvsLam") - # Normalise by good frames - GoodFrames = groupGet('IvsLam','samp','goodfrm') - print "run frames: ", GoodFrames - if (run=='0'): - RunNumber = '0' - else: - RunNumber = groupGet('IvsLam','samp','run_number') - #mantid['IvsLam'].getRun().getLogData("goodfrm").value - #Scale('IvsLam','IvsLam',GoodFrames**-1,'Multiply') - #IvsLam = mantid['IvsLam']*GoodFrames**-1 - if (trans==''): - #monitor2Eff('M') # This doesn't seem to work. - #heliumDetectorEff('DP') # point detector #Nor does this. - # Multidetector (Flood) TODO - print "No transmission file. Trying default exponential/polynomial correction..." - inst=groupGet('_DP','inst') - corrType=inst.getStringParameter('correction')[0] - if (corrType=='polynomial'): - pString=inst.getStringParameter('polystring') - print pString - if len(pString): - PolynomialCorrection(InputWorkspace='_DP',OutputWorkspace='IvsLam',Coefficients=pString[0],Operation='Divide') - else: - print "No polynomial coefficients in IDF. Using monitor spectrum with no corrections." - elif (corrType=='exponential'): - c0=inst.getNumberParameter('C0') - c1=inst.getNumberParameter('C1') - print "Exponential parameters: ", c0[0], c1[0] - if len(c0): - ExponentialCorrection(InputWorkspace='_DP',OutputWorkspace='IvsLam',C0=c0[0],C1=c1[0],Operation='Divide') - # normalise by monitor spectrum - # RebinToWorkspace(WorkspaceToRebin="_M",WorkspaceToMatch="_DP",OutputWorkspace="_M_M") - # CropWorkspace(InputWorkspace="M_M",OutputWorkspace="I0M",StartWorkspaceIndex=I0MonitorIndex) - # Divide(LHSWorkspace="DM",RHSWorkspace="I0M",OutputWorkspace="RM") - Divide(LHSWorkspace="IvsLam",RHSWorkspace="_I0P",OutputWorkspace="IvsLam") - else: # we have a transmission run - Integration(InputWorkspace="_I0P",OutputWorkspace="_monInt",RangeLower=str(intmin),RangeUpper=str(intmax)) - #scaling=1/mantid.getMatrixWorkspace('_monInt').dataY(0)[0] - Divide(LHSWorkspace="_DP",RHSWorkspace="_monInt",OutputWorkspace="IvsLam") - ##Divide(LHSWorkspace="_DP",RHSWorkspace="_I0P",OutputWorkspace="IvsLam") - names = mtd.getObjectNames() - if trans in names: - ##Divide(LHSWorkspace="_DP",RHSWorkspace="_I0P",OutputWorkspace="IvsLam") - RebinToWorkspace(WorkspaceToRebin=trans,WorkspaceToMatch="IvsLam",OutputWorkspace=trans) - ##IvsLam = mantid['IvsLam']*GoodFrames**-1 - Divide(LHSWorkspace="IvsLam",RHSWorkspace=trans,OutputWorkspace="IvsLam") - else: - transCorr(trans) - # Need to process the optional args to see what needs to be output and what division needs to be made - - # Convert to I vs Q - # check if detector in direct beam - if (theta == 0 or theta == ''): - if (theta == ''): - theta = 0 - print "given theta = ",theta - inst = groupGet('IvsLam','inst') - detLocation=inst.getComponentByName('point-detector').getPos() - sampleLocation=inst.getComponentByName('some-surface-holder').getPos() - detLocation=inst.getComponentByName('point-detector').getPos() - sample2detector=detLocation-sampleLocation # metres - source=inst.getSource() - beamPos = sampleLocation - source.getPos() - PI = 3.1415926535 - theta = inst.getComponentByName('point-detector').getTwoTheta(sampleLocation, beamPos)*180.0/PI/2.0 - print "Det location: ", detLocation, "Calculated theta = ",theta - if detLocation.getY() == 0: # detector is not in correct place - print "det at 0" - # Load corresponding NeXuS file - runno = '_' + str(run) - #templist.append(runno) - if type(run)==type(int()): - LoadNexus(Filename=run,OutputWorkspace=runno) - else: - LoadNexus(Filename=run.replace("raw","nxs",1),OutputWorkspace=runno) - # Get detector angle theta from NeXuS - theta = groupGet(runno,'samp','theta') - print 'Nexus file theta =', theta - IvsQ = l2q(mtd['IvsLam'], 'point-detector', theta) - else: - ConvertUnits(InputWorkspace='IvsLam',OutputWorkspace="IvsQ",Target="MomentumTransfer") - - else: - theta = float(theta) - IvsQ = l2q(mtd['IvsLam'], 'point-detector', theta) - - RenameWorkspace(InputWorkspace='IvsLam',OutputWorkspace=RunNumber+'_IvsLam') - RenameWorkspace(InputWorkspace='IvsQ',OutputWorkspace=RunNumber+'_IvsQ') - - # delete all temporary workspaces unless in debug mode (debug=1) + ''' - if debug == 0: - cleanup() - return mtd[RunNumber+'_IvsLam'], mtd[RunNumber+'_IvsQ'], theta - - - -def transCorr(transrun): - inst = groupGet("_W",'inst') - # Some beamline constants from IDF - intmin = inst.getNumberParameter('MonitorIntegralMin')[0] - intmax = inst.getNumberParameter('MonitorIntegralMax')[0] - if ',' in transrun: - slam = transrun.split(',')[0] - llam = transrun.split(',')[1] - print "Transmission runs: ", transrun - [I0MonitorIndex, MultiDetectorStart, nHist] = toLam(slam,'_'+slam) - CropWorkspace(InputWorkspace="_D_"+slam,OutputWorkspace="_D_"+slam,StartWorkspaceIndex=0,EndWorkspaceIndex=0) - [I0MonitorIndex, MultiDetectorStart, nHist] = toLam(llam,'_'+llam) - CropWorkspace(InputWorkspace="_D_"+llam,OutputWorkspace="_D_"+llam,StartWorkspaceIndex=0,EndWorkspaceIndex=0) - - RebinToWorkspace(WorkspaceToRebin="_M_"+llam,WorkspaceToMatch="_DP_"+llam,OutputWorkspace="_M_P_"+llam) - CropWorkspace(InputWorkspace="_M_P_"+llam,OutputWorkspace="_I0P_"+llam,StartWorkspaceIndex=I0MonitorIndex) - - #Normalise by monitor integral - inst = groupGet('_D_'+slam,'inst') - # Some beamline constants from IDF - intmin = inst.getNumberParameter('MonitorIntegralMin')[0] - intmax = inst.getNumberParameter('MonitorIntegralMax')[0] - Integration(InputWorkspace="_I0P_"+llam,OutputWorkspace="_monInt_TRANS",RangeLower=str(intmin),RangeUpper=str(intmax)) - Divide(LHSWorkspace="_DP_"+llam,RHSWorkspace="_monInt_TRANS",OutputWorkspace="_D_"+llam) - #scaling=1/mantid.getMatrixWorkspace('_monInt_TRANS').dataY(0)[0] - #Scale(InputWorkspace="_DP_"+llam,OutputWorkspace="_D_"+llam,Factor=scaling,Operation="Multiply") - - # same for short wavelength run slam: - RebinToWorkspace(WorkspaceToRebin="_M_"+slam,WorkspaceToMatch="_DP_"+slam,OutputWorkspace="_M_P_"+slam) - CropWorkspace(InputWorkspace="_M_P_"+slam,OutputWorkspace="_I0P_"+slam,StartWorkspaceIndex=I0MonitorIndex) + run_ws = ConvertToWavelength.to_workspace(run) + idf_defaults = get_defaults(run_ws) + to_lam = ConvertToWavelength(run_ws) + nHist = run_ws.getNumberHistograms() + + I0MonitorIndex = idf_defaults['I0MonitorIndex'] + MultiDetectorStart = idf_defaults['MultiDetectorStart'] + lambda_min = idf_defaults['LambdaMin'] + lambda_max = idf_defaults['LambdaMax'] + detector_index_ranges = (idf_defaults['PointDetectorStart'], idf_defaults['PointDetectorStop']) + background_min = idf_defaults['MonitorBackgroundMin'] + background_max = idf_defaults['MonitorBackgroundMax'] + intmin = idf_defaults['MonitorIntegralMin'] + intmax = idf_defaults['MonitorIntegralMax'] + + _monitor_ws, _detector_ws = to_lam.convert(wavelength_min=lambda_min, wavelength_max=lambda_max, detector_workspace_indexes=detector_index_ranges, monitor_workspace_index=I0MonitorIndex, correct_monitor=True, bg_min=background_min, bg_max=background_max ) + + inst = run_ws.getInstrument() + # Some beamline constants from IDF + + print I0MonitorIndex + print nHist + if (nHist > 5 and not(pointdet)): + # Proccess Multi-Detector; assume MD goes to the end: + # if roi or db are given in the function then sum over the apropriate channels + print "This is a multidetector run." + try: + _I0M = RebinToWorkspace(WorkspaceToRebin=_monitor_ws,WorkspaceToMatch=_detector_ws) + IvsLam = _detector_ws / _IOM + if (roi != [0,0]) : + ReflectedBeam = SumSpectra(InputWorkspace=IvsLam, StartWorkspaceIndex=roi[0], EndWorkspaceIndex=roi[1]) + if (db != [0,0]) : + DirectBeam = SumSpectra(InputWorkspace=_detector_ws, StartWorkspaceIndex=db[0], EndWorkspaceIndex=db[1]) + except SystemExit: + print "Point-Detector only run." + RunNumber = groupGet(IvsLam.getName(),'samp','run_number') + if (theta): + IvsQ = l2q(ReflectedBeam, 'linear-detector', theta) # TODO: possible to get here and an invalid state if roi == [0,0] see above. + else: + IvsQ = ConvertUnits(InputWorkspace=ReflectedBeam, Target="MomentumTransfer") + + # Single Detector processing------------------------------------------------------------- + else: + print "This is a Point-Detector run." + # handle transmission runs + # process the point detector reflectivity + _I0P = RebinToWorkspace(WorkspaceToRebin=_monitor_ws,WorkspaceToMatch=_detector_ws) + IvsLam = Scale(InputWorkspace=_detector_ws,Factor=1) + # Normalise by good frames + GoodFrames = groupGet(IvsLam.getName(),'samp','goodfrm') + print "run frames: ", GoodFrames + if (run=='0'): + RunNumber = '0' + else: + RunNumber = groupGet(IvsLam.getName(),'samp','run_number') + if (trans==''): + print "No transmission file. Trying default exponential/polynomial correction..." + inst=groupGet(_detector_ws.getName(),'inst') + corrType=inst.getStringParameter('correction')[0] + if (corrType=='polynomial'): + pString=inst.getStringParameter('polystring') + print pString + if len(pString): + IvsLam = PolynomialCorrection(InputWorkspace=_detector_ws,Coefficients=pString[0],Operation='Divide') + else: + print "No polynomial coefficients in IDF. Using monitor spectrum with no corrections." + elif (corrType=='exponential'): + c0=inst.getNumberParameter('C0') + c1=inst.getNumberParameter('C1') + print "Exponential parameters: ", c0[0], c1[0] + if len(c0): + IvsLam = ExponentialCorrection(InputWorkspace=_detector_ws,C0=c0[0],C1=c1[0],Operation='Divide') + IvsLam = Divide(LHSWorkspace=IvsLam, RHSWorkspace=_I0P) + else: # we have a transmission run + _monInt = Integration(InputWorkspace=_I0P,RangeLower=intmin,RangeUpper=intmax) + IvsLam = Divide(LHSWorkspace=_detector_ws,RHSWorkspace=_monInt) + names = mtd.getObjectNames() + if trans in names: + trans = RebinToWorkspace(WorkspaceToRebin=trans,WorkspaceToMatch=IvsLam,OutputWorkspace=trans) + IvsLam = Divide(LHSWorkspace=IvsLam,RHSWorkspace=trans,OutputWorkspace="IvsLam") # TODO: Hardcoded names are bad + else: + IvsLam = transCorr(trans, IvsLam) + print type(IvsLam) + RenameWorkspace(InputWorkspace=IvsLam, OutputWorkspace="IvsLam") # TODO: Hardcoded names are bad + + + # Convert to I vs Q + # check if detector in direct beam + if (theta == 0 or theta == ''): + if (theta == ''): + theta = 0 + print "given theta = ",theta + inst = groupGet('IvsLam','inst') + detLocation=inst.getComponentByName('point-detector').getPos() + sampleLocation=inst.getComponentByName('some-surface-holder').getPos() + detLocation=inst.getComponentByName('point-detector').getPos() + sample2detector=detLocation-sampleLocation # metres + source=inst.getSource() + beamPos = sampleLocation - source.getPos() + PI = 3.1415926535 + theta = inst.getComponentByName('point-detector').getTwoTheta(sampleLocation, beamPos)*180.0/PI/2.0 + print "Det location: ", detLocation, "Calculated theta = ",theta + if detLocation.getY() == 0: # detector is not in correct place + print "det at 0" + # Load corresponding NeXuS file + runno = '_' + str(run) + #templist.append(runno) + if type(run)==type(int()): + LoadNexus(Filename=run,OutputWorkspace=runno) + else: + LoadNexus(Filename=run.replace("raw","nxs",1),OutputWorkspace=runno) + # Get detector angle theta from NeXuS + theta = groupGet(runno,'samp','theta') + print 'Nexus file theta =', theta + IvsQ = l2q(mtd['IvsLam'], 'point-detector', theta) + else: + ConvertUnits(InputWorkspace='IvsLam',OutputWorkspace="IvsQ",Target="MomentumTransfer") + + else: + theta = float(theta) + IvsQ = l2q(mtd['IvsLam'], 'point-detector', theta) + + RenameWorkspace(InputWorkspace='IvsLam',OutputWorkspace=RunNumber+'_IvsLam') + RenameWorkspace(InputWorkspace='IvsQ',OutputWorkspace=RunNumber+'_IvsQ') + + # delete all temporary workspaces unless in debug mode (debug=1) + + if debug == 0: + pass + cleanup() + return mtd[RunNumber+'_IvsLam'], mtd[RunNumber+'_IvsQ'], theta - #Normalise by monitor integral - inst = groupGet('_D_'+llam,'inst') - # Some beamline constants from IDF - intmin = inst.getNumberParameter('MonitorIntegralMin')[0] - intmax = inst.getNumberParameter('MonitorIntegralMax')[0] - Integration(InputWorkspace="_I0P_"+slam,OutputWorkspace="_monInt_TRANS",RangeLower=str(intmin),RangeUpper=str(intmax)) - #scaling=1/mantid.getMatrixWorkspace('_monInt_TRANS').dataY(0)[0] - Divide(LHSWorkspace="_DP_"+slam,RHSWorkspace="_monInt_TRANS",OutputWorkspace="_D_"+slam) - #Scale(InputWorkspace="_DP_"+slam,OutputWorkspace="_D_"+slam,Factor=scaling,Operation="Multiply") - - #Divide(LHSWorkspace="_DP_"+slam,RHSWorkspace="_I0P_"+slam,OutputWorkspace="_D_"+slam) - [transr, sf] = combine2("_D_"+slam,"_D_"+llam,"_DP_TRANS",10.0,12.0,1.5,17.0,0.02,scalehigh=1) - #[wlam, wq, th] = quick(runno,angle,trans='_transcomb') - else: - [I0MonitorIndex, MultiDetectorStart, nHist] = toLam(transrun,'_TRANS') - RebinToWorkspace(WorkspaceToRebin="_M_TRANS",WorkspaceToMatch="_DP_TRANS",OutputWorkspace="_M_P_TRANS") - CropWorkspace(InputWorkspace="_M_P_TRANS",OutputWorkspace="_I0P_TRANS",StartWorkspaceIndex=I0MonitorIndex) - #Normalise by monitor integral - Integration(InputWorkspace="_I0P_TRANS",OutputWorkspace="_monInt_TRANS",RangeLower=str(intmin),RangeUpper=str(intmax)) - Divide(LHSWorkspace="_DP_TRANS",RHSWorkspace="_monInt_TRANS",OutputWorkspace="_DP_TRANS") - #scaling=1/mantid.getMatrixWorkspace('_monInt_TRANS').dataY(0)[0] - #print "SCALING:",scaling - #Scale(InputWorkspace="_I0P_TRANS",OutputWorkspace=str(transrun)+"_IvsLam_TRANS",Factor=scaling,Operation="Multiply") - #Scale(InputWorkspace="_DP_TRANS",OutputWorkspace="_DP_TRANS",Factor=scaling,Operation="Multiply") - - #got sometimes very slight binning diferences, so do this again: - RebinToWorkspace(WorkspaceToRebin='_DP_TRANS',WorkspaceToMatch="IvsLam",OutputWorkspace=str(transrun)+'_IvsLam_TRANS') - if isinstance(mtd["_DP_TRANS"], WorkspaceGroup): - Divide(LHSWorkspace="IvsLam",RHSWorkspace=str(transrun)+"_IvsLam_TRANS_1",OutputWorkspace="IvsLam") - else: - Divide(LHSWorkspace="IvsLam",RHSWorkspace=str(transrun)+"_IvsLam_TRANS",OutputWorkspace="IvsLam") +def transCorr(transrun, i_vs_lam): + + run_ws = ConvertToWavelength.to_workspace(i_vs_lam) + idf_defaults = get_defaults(run_ws) + I0MonitorIndex = idf_defaults['I0MonitorIndex'] + MultiDetectorStart = idf_defaults['MultiDetectorStart'] + lambda_min = idf_defaults['LambdaMin'] + lambda_max = idf_defaults['LambdaMax'] + background_min = idf_defaults['MonitorBackgroundMin'] + background_max = idf_defaults['MonitorBackgroundMax'] + intmin = idf_defaults['MonitorIntegralMin'] + intmax = idf_defaults['MonitorIntegralMax'] + detector_index_ranges = (idf_defaults['PointDetectorStart'], idf_defaults['PointDetectorStop']) + + + transWS = None + if ',' in transrun: + slam = transrun.split(',')[0] + llam = transrun.split(',')[1] + print "Transmission runs: ", transrun + + to_lam = ConvertToWavelength(slam) + _monitor_ws_slam, _detector_ws_slam = to_lam.convert(wavelength_min=lambda_min, wavelength_max=lambda_max, detector_workspace_indexes=detector_index_ranges, monitor_workspace_index=I0MonitorIndex, correct_monitor=True, bg_min=background_min, bg_max=background_max ) + + _i0p_slam = RebinToWorkspace(WorkspaceToRebin=_monitor_ws_slam, WorkspaceToMatch=_detector_ws_slam) + _mon_int_trans = Integration(InputWorkspace=_i0p_slam, RangeLower=intmin, RangeUpper=intmax) + _detector_ws_slam = Divide(LHSWorkspace=_detector_ws_slam, RHSWorkspace=_mon_int_trans) + + to_lam = ConvertToWavelength(llam) + _monitor_ws_llam, _detector_ws_llam = to_lam.convert(wavelength_min=lambda_min, wavelength_max=lambda_max, detector_workspace_indexes=detector_index_ranges, monitor_workspace_index=I0MonitorIndex, correct_monitor=True, bg_min=background_min, bg_max=background_max ) + + _i0p_llam = RebinToWorkspace(WorkspaceToRebin=_monitor_ws_llam, WorkspaceToMatch=_detector_ws_llam) + _mon_int_trans = Integration(InputWorkspace=_i0p_llam, RangeLower=intmin,RangeUpper=intmax) + _detector_ws_llam = Divide(LHSWorkspace=_detector_ws_llam, RHSWorkspace=_mon_int_trans) + + # TODO: HARDCODED STITCHING VALUES!!!!! + _transWS, outputScaling = Stitch1D(LHSWorkspace=_detector_ws_slam, RHSWorkspace=_detector_ws_llam, StartOverlap=10, EndOverlap=12, Params="%f,%f,%f" % (1.5, 0.02, 17)) + + else: + + to_lam = ConvertToWavelength(transrun) + _monitor_ws_trans, _detector_ws_trans = to_lam.convert(wavelength_min=lambda_min, wavelength_max=lambda_max, detector_workspace_indexes=detector_index_ranges, monitor_workspace_index=I0MonitorIndex, correct_monitor=True, bg_min=background_min, bg_max=background_max ) + _i0p_trans = RebinToWorkspace(WorkspaceToRebin=_monitor_ws_trans, WorkspaceToMatch=_detector_ws_trans) + + _mon_int_trans = Integration( InputWorkspace=_i0p_trans, RangeLower=intmin, RangeUpper=intmax ) + _transWS = Divide( LHSWorkspace=_detector_ws_trans, RHSWorkspace=_mon_int_trans ) + + + #got sometimes very slight binning diferences, so do this again: + _i_vs_lam_trans = RebinToWorkspace(WorkspaceToRebin=_transWS, WorkspaceToMatch=i_vs_lam) + # Normalise by transmission run. + _i_vs_lam_corrected = i_vs_lam / _i_vs_lam_trans + + return _i_vs_lam_corrected def cleanup(): - names = mtd.getObjectNames() - for name in names: - if re.search("^_", name): - DeleteWorkspace(name) - -def coAdd(run,name): - names = mtd.getObjectNames() - wTof = "_W" + name # main workspace in time-of-flight - if run in names: - RenameWorkspace(InputWorkspace=run,OutputWorkspace=wTof) - else: - - currentInstrument=config['default.instrument'] - runlist = [] - l1 = run.split(',') - for subs in l1: - l2 = subs.split(':') - for l3 in l2: - runlist.append(l3) - print "Adding: ", runlist - currentInstrument=currentInstrument.upper() - if (runlist[0]=='0'): #DAE/current run - StartLiveData(Instrument=currentInstrument,UpdateEvery='0',Outputworkspace='_sum') - #LoadLiveData(currentInstrument,OutputWorkspace='_sum') - #LoadDAE(DAEname='ndx'+mantid.settings['default.instrument'],OutputWorkspace='_sum') - else: - Load(Filename=runlist[0],OutputWorkspace='_sum')#,LoadLogFiles="1") - - for i in range(len(runlist)-1): - if (runlist[i+1]=='0'): #DAE/current run - StartLiveData(Instrument=currentInstrument,UpdateEvery='0',Outputworkspace='_w2') - #LoadLiveData(currentInstrument,OutputWorkspace='_w2') - #LoadDAE(DAEname='ndx'+mantid.settings['default.instrument'],OutputWorkspace='_w2') - else: - Load(Filename=runlist[i+1],OutputWorkspace='_w2')#,LoadLogFiles="1") - - Plus(LHSWorkspace='_sum',RHSWorkspace='_w2',OutputWorkspace='_sum') - - RenameWorkspace(InputWorkspace='_sum',OutputWorkspace=wTof) - - - -def toLam(run, name, pointdet=1): - ''' - toLam splits a given run into monitor and detector spectra and - converts these to wavelength - ''' - # some naming for convenience - wTof = "_W" + name # main workspace in time-of-flight - monInLam = "_M" + name # monitor spectra vs. wavelength - detInLam = "_D" + name # detector spectra vs. wavelength - pDet = "_DP" + name # point-detector only vs. wavelength - mDet = "_DM" + name # multi-detector only vs. wavelength - - - # add multiple workspaces, if given - coAdd(run,name) - - inst = groupGet(wTof,'inst') - # Some beamline constants from IDF - - bgmin = inst.getNumberParameter('MonitorBackgroundMin')[0] - bgmax = inst.getNumberParameter('MonitorBackgroundMax')[0] - MonitorBackground = [bgmin,bgmax] - intmin = inst.getNumberParameter('MonitorIntegralMin')[0] - intmax = inst.getNumberParameter('MonitorIntegralMax')[0] - MonitorsToCorrect = [int(inst.getNumberParameter('MonitorsToCorrect')[0])] - # Note: Since we are removing the monitors in the load raw command they are not counted here. - PointDetectorStart = int(inst.getNumberParameter('PointDetectorStart')[0]) - PointDetectorStop = int(inst.getNumberParameter('PointDetectorStop')[0]) - MultiDetectorStart = int(inst.getNumberParameter('MultiDetectorStart')[0]) - I0MonitorIndex = int(inst.getNumberParameter('I0MonitorIndex')[0]) - - # Get usable wavelength range - LambdaMin = float(inst.getNumberParameter('LambdaMin')[0]) - LambdaMax = float(inst.getNumberParameter('LambdaMax')[0]) - - # Convert spectra from TOF to wavelength - ConvertUnits(InputWorkspace=wTof,OutputWorkspace=wTof+"_lam",Target="Wavelength", AlignBins='1') - # Separate detector an monitor spectra manually - CropWorkspace(InputWorkspace=wTof+"_lam",OutputWorkspace=monInLam,StartWorkspaceIndex='0',EndWorkspaceIndex=PointDetectorStart-1) - CropWorkspace(InputWorkspace=wTof+"_lam",OutputWorkspace=detInLam,XMin=LambdaMin,XMax=LambdaMax,StartWorkspaceIndex=PointDetectorStart) - # Subtract flat background from fit in range given from Instrument Def/Par File - CalculateFlatBackground(InputWorkspace=monInLam,OutputWorkspace=monInLam,WorkspaceIndexList=MonitorsToCorrect,StartX=MonitorBackground[0],EndX=MonitorBackground[1]) - - # Is it a multidetector run? - nHist = groupGet(wTof+"_lam",'wksp','') - if (nHist<6 or (nHist>5 and pointdet)): - CropWorkspace(InputWorkspace=wTof+"_lam",OutputWorkspace=pDet,XMin=LambdaMin,XMax=LambdaMax,StartWorkspaceIndex=PointDetectorStart,EndWorkspaceIndex=PointDetectorStop) - else: - CropWorkspace(InputWorkspace=wTof+"_lam",OutputWorkspace=mDet,XMin=LambdaMin,XMax=LambdaMax,StartWorkspaceIndex=MultiDetectorStart) - - - - - return I0MonitorIndex, MultiDetectorStart, nHist - + names = mtd.getObjectNames() + for name in names: + if re.search("^_", name): + DeleteWorkspace(name) + +def get_defaults(run_ws): + ''' + Temporary helper function. Aid refactoring by removing need to specifically ask things of parameter files. + ''' + defaults = dict() + if isinstance(run_ws, WorkspaceGroup): + instrument = run_ws[0].getInstrument() + else: + instrument = run_ws.getInstrument() + defaults['LambdaMin'] = float( instrument.getNumberParameter('LambdaMin')[0] ) + defaults['LambdaMax'] = float( instrument.getNumberParameter('LambdaMax')[0] ) + defaults['MonitorBackgroundMin'] = float( instrument.getNumberParameter('MonitorBackgroundMin')[0] ) + defaults['MonitorBackgroundMax'] = float( instrument.getNumberParameter('MonitorBackgroundMax')[0] ) + defaults['MonitorIntegralMin'] = float( instrument.getNumberParameter('MonitorIntegralMin')[0] ) + defaults['MonitorIntegralMax'] = float( instrument.getNumberParameter('MonitorIntegralMax')[0] ) + defaults['MonitorsToCorrect'] = int( instrument.getNumberParameter('MonitorsToCorrect')[0] ) + defaults['PointDetectorStart'] = int( instrument.getNumberParameter('PointDetectorStart')[0] ) + defaults['PointDetectorStop'] = int( instrument.getNumberParameter('PointDetectorStop')[0] ) + defaults['MultiDetectorStart'] = int( instrument.getNumberParameter('MultiDetectorStart')[0] ) + defaults['I0MonitorIndex'] = int( instrument.getNumberParameter('I0MonitorIndex')[0] ) + return defaults + def groupGet(wksp,whattoget,field=''): - ''' - returns information about instrument or sample details for a given workspace wksp, - also if the workspace is a group (info from first group element) - ''' - if (whattoget == 'inst'): - if isinstance(mtd[wksp], WorkspaceGroup): - return mtd[wksp+'_1'].getInstrument() - else: - return mtd[wksp].getInstrument() - - elif (whattoget == 'samp' and field != ''): - if isinstance(mtd[wksp], WorkspaceGroup): - try: - log = mtd[wksp + '_1'].getRun().getLogData(field).value - if (type(log) is int or type(log) is str): - res=log - else: - res = log[len(log)-1] - except RuntimeError: - res = 0 - print "Block "+field+" not found." - else: - try: - log = mtd[wksp].getRun().getLogData(field).value - if (type(log) is int or type(log) is str): - res=log - else: - res = log[len(log)-1] - except RuntimeError: - res = 0 - print "Block "+field+" not found." - return res - elif (whattoget == 'wksp'): - if isinstance(mtd[wksp], WorkspaceGroup): - return mtd[wksp+'_1'].getNumberHistograms() - else: - return mtd[wksp].getNumberHistograms() + ''' + returns information about instrument or sample details for a given workspace wksp, + also if the workspace is a group (info from first group element) + ''' + if (whattoget == 'inst'): + if isinstance(mtd[wksp], WorkspaceGroup): + return mtd[wksp+'_1'].getInstrument() + else: + return mtd[wksp].getInstrument() + + elif (whattoget == 'samp' and field != ''): + if isinstance(mtd[wksp], WorkspaceGroup): + try: + log = mtd[wksp + '_1'].getRun().getLogData(field).value + if (type(log) is int or type(log) is str): + res=log + else: + res = log[len(log)-1] + except RuntimeError: + res = 0 + print "Block "+field+" not found." + else: + try: + log = mtd[wksp].getRun().getLogData(field).value + if (type(log) is int or type(log) is str): + res=log + else: + res = log[len(log)-1] + except RuntimeError: + res = 0 + print "Block "+field+" not found." + return res + elif (whattoget == 'wksp'): + if isinstance(mtd[wksp], WorkspaceGroup): + return mtd[wksp+'_1'].getNumberHistograms() + else: + return mtd[wksp].getNumberHistograms() """ ===================== Testing stuff ===================== - - ToDo: - 1) Make the test below more rigorous. - 2) Each test should either print out Success or Failure - or describe in words what the tester should see on the screen. - 3) Each function should be accompanied by a test function - 4) all test functions to be run in a give file should be called in doAllTests(). - The reason for this is that the normal python if __name__ == '__main__': - testing location is heavily used during script development and debugging. + + ToDo: + 1) Make the test below more rigorous. + 2) Each test should either print out Success or Failure + or describe in words what the tester should see on the screen. + 3) Each function should be accompanied by a test function + 4) all test functions to be run in a give file should be called in doAllTests(). + The reason for this is that the normal python if __name__ == '__main__': + testing location is heavily used during script development and debugging. """ - + def _testQuick(): - config['default.instrument'] = "SURF" - [w1lam,w1q,th] = quick(94511,theta=0.25,trans='94504') - [w2lam,w2q,th] = quick(94512,theta=0.65,trans='94504') - [w3lam,w3q,th] = quick(94513,theta=1.5,trans='94504') - g1=plotSpectrum("94511_IvsQ",0) - g2=plotSpectrum("94512_IvsQ",0) - g3=plotSpectrum("94513_IvsQ",0) - - return True + config['default.instrument'] = "SURF" + [w1lam,w1q,th] = quick(94511,theta=0.25,trans='94504') + [w2lam,w2q,th] = quick(94512,theta=0.65,trans='94504') + [w3lam,w3q,th] = quick(94513,theta=1.5,trans='94504') + g1=plotSpectrum("94511_IvsQ",0) + g2=plotSpectrum("94512_IvsQ",0) + g3=plotSpectrum("94513_IvsQ",0) + + return True def _doAllTests(): - _testQuick() - return True + _testQuick() + return True if __name__ == '__main__': - ''' This is the debugging and testing area of the file. The code below is run when ever the - script is called directly from a shell command line or the execute all menu option in mantid. - ''' - #Debugging = True # Turn the debugging on and the testing code off - Debugging = False # Turn the debugging off and the testing on - - if Debugging == False: - _doAllTests() - else: #Debugging code goes below - rr = quick("N:/SRF93080.raw") - + ''' This is the debugging and testing area of the file. The code below is run when ever the + script is called directly from a shell command line or the execute all menu option in mantid. + ''' + #Debugging = True # Turn the debugging on and the testing code off + Debugging = False # Turn the debugging off and the testing on + + if Debugging == False: + _doAllTests() + else: #Debugging code goes below + rr = quick("N:/SRF93080.raw") + # x=quick(95266) diff --git a/Code/Mantid/scripts/test/ConvertToWavelengthTest.py b/Code/Mantid/scripts/test/ConvertToWavelengthTest.py new file mode 100644 index 000000000000..5334d0a58835 --- /dev/null +++ b/Code/Mantid/scripts/test/ConvertToWavelengthTest.py @@ -0,0 +1,120 @@ +import unittest +from mantid.simpleapi import * +from isis_reflectometry.convert_to_wavelength import ConvertToWavelength + +class ConvertToWavelengthTest(unittest.TestCase): + """ + Test the convert to wavelength type. + """ + def test_construction_from_single_ws(self): + ws = CreateWorkspace(DataY=[1,2,3], DataX=[1,2,3]) + converter = ConvertToWavelength(ws) + self.assertTrue(converter != None, "Should have been able to make a valid converter from a single workspace") + DeleteWorkspace(ws) + + def test_construction_from_single_ws_name(self): + ws = CreateWorkspace(DataY=[1,2,3], DataX=[1,2,3]) + + converter = ConvertToWavelength(ws.getName()) + self.assertTrue(converter != None, "Should have been able to make a valid converter from a single workspace name") + DeleteWorkspace(ws) + + def test_construction_from_many_workspaces(self): + ws1 = CreateWorkspace(DataY=[1,2,3], DataX=[1,2,3]) + ws2 = CreateWorkspace(DataY=[1,2,3], DataX=[1,2,3]) + converter = ConvertToWavelength([ws1, ws2]) + self.assertTrue(converter != None, "Should have been able to make a valid converter from many workspace objects") + DeleteWorkspace(ws1) + DeleteWorkspace(ws2) + + def test_construction_from_many_workspace_names(self): + ws1 = CreateWorkspace(DataY=[1,2,3], DataX=[1,2,3]) + ws2 = CreateWorkspace(DataY=[1,2,3], DataX=[1,2,3]) + converter = ConvertToWavelength([ws1.getName(), ws2.getName()]) + self.assertTrue(converter != None, "Should have been able to make a valid converter from many workspace objects") + DeleteWorkspace(ws1) + DeleteWorkspace(ws2) + + def test_sum_workspaces(self): + ws1 = CreateWorkspace(DataY=[1,2,3], DataX=[1,2,3]) + ws2 = CloneWorkspace(ws1) + ws3 = CloneWorkspace(ws1) + sum = ConvertToWavelength.sum_workspaces([ws1, ws2, ws3]) + self.assertEqual(set([3,6,9]), set(sum.readY(0)), "Fail to sum workspaces correctly") + DeleteWorkspace(ws1) + DeleteWorkspace(ws2) + DeleteWorkspace(ws3) + DeleteWorkspace(sum) + + def test_conversion_throws_with_min_wavelength_greater_or_equal_to_max_wavelength(self): + ws = CreateWorkspace(DataY=[1,2,3], DataX=[1,2,3]) + converter = ConvertToWavelength(ws) + self.assertRaises(ValueError, converter.convert, 1.0, 0.0, (), 0) + self.assertRaises(ValueError, converter.convert, 1.0, 1.0, (), 0) + DeleteWorkspace(ws) + + def test_conversion_throws_with_some_flat_background_params_but_not_all(self): + ws = CreateWorkspace(DataY=[1,2,3], DataX=[1,2,3]) + converter = ConvertToWavelength(ws) + self.assertRaises(ValueError, converter.convert, 0.0, 1.0, (), 0, True) + DeleteWorkspace(ws) + + def test_conversion_throws_with_min_background_greater_than_or_equal_to_max_background(self): + ws = CreateWorkspace(DataY=[1,2,3], DataX=[1,2,3]) + converter = ConvertToWavelength(ws) + self.assertRaises(ValueError, converter.convert, 0.0, 1.0, (), 0, True, 1.0, 1.0) + DeleteWorkspace(ws) + + + def test_crop_range(self): + original_ws = Load(Filename='INTER00013460') + + # Crop out one spectra + temp_ws = ConvertToWavelength.crop_range(original_ws, (0, original_ws.getNumberHistograms()-2)) + self.assertEqual(original_ws.getNumberHistograms()-1, temp_ws.getNumberHistograms()) + + # Crop out all but 2 spectra from start and end. + temp_ws = ConvertToWavelength.crop_range(original_ws, ( (0, 1), (3, 4) ) ) + self.assertEqual(2, temp_ws.getNumberHistograms()) + + # Crop out all but 2 spectra from start and end. Exactly the same as above, but slightly different tuple syntax + temp_ws = ConvertToWavelength.crop_range(original_ws, ( ( (0, 1), (3, 4) ) )) + self.assertEqual(2, temp_ws.getNumberHistograms()) + + # Crop out all but 1 spectra + temp_ws = ConvertToWavelength.crop_range(original_ws, ( (1, 3) ) ) + self.assertEqual(3, temp_ws.getNumberHistograms()) + # First and last dectors are cropped off, so indexes go 2-4 rather than 1-5 + self.assertEqual(2, temp_ws.getDetector(0).getID()) + self.assertEqual(3, temp_ws.getDetector(1).getID()) + self.assertEqual(4, temp_ws.getDetector(2).getID()) + + # Test resilience to junk inputs + self.assertRaises(ValueError, ConvertToWavelength.crop_range, original_ws, 'a') + self.assertRaises(ValueError, ConvertToWavelength.crop_range, original_ws, (1,2,3)) + + @classmethod + def cropped_x_range(cls, ws, index): + det_ws_x = ws.readX(index) + mask = ws.readY(index) != 0 # CropWorkspace will only zero out y values! so we need to translate those to an x range + cropped_x = det_ws_x[mask] + return cropped_x[0], cropped_x[-1] + + def test_convert(self): + ws = Load(Filename='INTER00013460') + converter = ConvertToWavelength(ws) + + monitor_ws, detector_ws = converter.convert(wavelength_min=0.0, wavelength_max=10.0, detector_workspace_indexes = (2,4), monitor_workspace_index=0, correct_monitor=True, bg_min=2.0, bg_max=8.0) + + print detector_ws.getNumberHistograms() + self.assertEqual(1, monitor_ws.getNumberHistograms(), "Wrong number of spectra in monitor workspace") + self.assertEqual(3, detector_ws.getNumberHistograms(), "Wrong number of spectra in detector workspace") + self.assertEqual("Wavelength", detector_ws.getAxis(0).getUnit().unitID()) + self.assertEqual("Wavelength", monitor_ws.getAxis(0).getUnit().unitID()) + x_min, x_max = ConvertToWavelengthTest.cropped_x_range(detector_ws, 0) + + self.assertTrue(x_min >= 0.0) + self.assertTrue(x_max <= 10.0) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/Code/Mantid/scripts/test/ReflectometryQuickAuxiliaryTest.py b/Code/Mantid/scripts/test/ReflectometryQuickAuxiliaryTest.py index 24e3695563eb..7eb579787cea 100644 --- a/Code/Mantid/scripts/test/ReflectometryQuickAuxiliaryTest.py +++ b/Code/Mantid/scripts/test/ReflectometryQuickAuxiliaryTest.py @@ -30,44 +30,6 @@ def test_cleanup(self): DeleteWorkspace(tokeep) - def test_coAdd_ws_in_ADS(self): - inWS = CreateSingleValuedWorkspace(DataValue=1, ErrorValue=1) - quick.coAdd('inWS', 'ProvidedName') - outWS = mtd['_WProvidedName'] - result = CheckWorkspacesMatch(Workspace1=inWS, Workspace2=outWS) - self.assertEquals("Success!", result) - DeleteWorkspace(outWS) - - def test_coAdd_run_list(self): - originalInstrument = config.getInstrument() - try: - # We have multiple runs from some MUSR files in AutoTest, lets use those. - tempInstrument = "MUSR" - config['default.instrument'] = tempInstrument - runlist = '15189, 15190' - - # Run coAdd - quick.coAdd(runlist, 'ProvidedName') - - # Get the output workspace and do some quick sanity checks - outWS = mtd['_WProvidedName'] - self.assertEquals(outWS[0].getInstrument().getName(), tempInstrument) - - # Perform the addition of the two files manually - a = LoadMuonNexus(Filename='15189') - b = LoadMuonNexus(Filename='15190') - c = Plus(LHSWorkspace=a[0], RHSWorkspace=b[0]) - - #Check the expected calculated result against coAdd - result = CheckWorkspacesMatch(Workspace1=c, Workspace2=outWS) - self.assertEquals("Success!", result) - finally: - config['default.instrument'] = originalInstrument.name() - DeleteWorkspace(a[0]) - DeleteWorkspace(b[0]) - DeleteWorkspace(c) - DeleteWorkspace(outWS) - def test_groupGet_instrument(self): expectedInstrument = "POLREF" diff --git a/Code/Mantid/scripts/test/ReflectometryQuickToLamTest.py b/Code/Mantid/scripts/test/ReflectometryQuickToLamTest.py deleted file mode 100644 index 8d82e502f7ab..000000000000 --- a/Code/Mantid/scripts/test/ReflectometryQuickToLamTest.py +++ /dev/null @@ -1,89 +0,0 @@ -import unittest -import numpy -from mantid.simpleapi import * -from mantid.api import * -from isis_reflectometry import quick - -class ReflectometryQuickToLamTest(unittest.TestCase): - - __wsName = None - - def setUp(self): - self.__wsName = "TestWorkspace" - LoadISISNexus(Filename='INTER00013460', OutputWorkspace=self.__wsName) - - def tearDown(self): - if self.__wsName in mtd.getObjectNames(): - DeleteWorkspace(mtd[self.__wsName]) - - def test_basic_type_checks(self): - firstWSName = mtd[self.__wsName].name() - - # Run Quick - quick.toLam(firstWSName, firstWSName) - - object_names = mtd.getObjectNames() - self.assertTrue('_W'+self.__wsName in object_names) - self.assertTrue('_M'+self.__wsName in object_names) - - ''' Internal renaming of inputs. Quick will need refactoring to stop this sort of thing, but - at least it get's documented here. - ''' - detectorWSInTOF = mtd['_W' + self.__wsName] - monitiorWSInTOF = mtd['_M' + self.__wsName] - # Fetch output workspace - detectorWSInLam = mtd['_W' + self.__wsName + '_lam'] - # Check workspace type - self.assertTrue(isinstance(detectorWSInLam, MatrixWorkspace)) - # Check units have switched to wavelength (hence lam - lambda in name) - self.assertEquals("Wavelength", detectorWSInLam.getAxis(0).getUnit().caption()) - - - def test_check_lambda_range(self): - - # Expected Min and Max x wavelength values on output. - inst = mtd[self.__wsName].getInstrument() - expectedLambdaMin = inst.getNumberParameter('LambdaMin')[0] - expectedLambdaMax = inst.getNumberParameter('LambdaMax')[0] - - firstWSName = mtd[self.__wsName].name() - - # Run Quick - quick.toLam(firstWSName, firstWSName) - - # Get output workspace - detectorWSInLam = mtd['_D' + self.__wsName] - - # Check that output workspace has been cropped. - x = detectorWSInLam.readX(0) - self.assertAlmostEquals(expectedLambdaMin, x[0], 0) - self.assertAlmostEquals(expectedLambdaMax, x[-1], 0) - - def test_workspace_splitting_monitor_detector(self): - # Expected Min and Max x wavelength values on output. - inst = mtd[self.__wsName].getInstrument() - originalNHisto = mtd[self.__wsName].getNumberHistograms() - - pointDetectorStartIndex = inst.getNumberParameter('PointDetectorStart')[0] - pointDetectorStopIndex = inst.getNumberParameter('PointDetectorStop')[0] - multiDetectorStartIndex = inst.getNumberParameter('MultiDetectorStart')[0] - - histoRangeMonitor = pointDetectorStartIndex # zero to pointDetectorStartIndex - histoRangePointDetector = originalNHisto - pointDetectorStartIndex - histoRangeMultiDetector = originalNHisto - multiDetectorStartIndex - - firstWSName = mtd[self.__wsName].name() - - # Run Quick - quick.toLam(firstWSName, firstWSName) - - # Get output workspace - pointDetectorWSInLam = mtd['_D' + self.__wsName] - monitorWSInLam = mtd['_M' + self.__wsName] - - # Check histogram ranges - self.assertEquals(histoRangeMonitor, monitorWSInLam.getNumberHistograms()) - self.assertEquals(histoRangePointDetector, pointDetectorWSInLam.getNumberHistograms()) - -if __name__ == '__main__': - unittest.main() \ No newline at end of file