Ticket #8643: quick_Max.py

File quick_Max.py, 28.8 KB (added by Owen Arnold, 7 years ago)
Line 
1''' SVN Info:   The variables below will only get subsituted at svn checkout if
2                the repository is configured for variable subsitution.
3
4        $Id$
5        $HeadURL$
6|=============================================================================|=======|
71                                                                            80   <tab>
8'''
9#these need to be moved into one NR folder or so
10#from ReflectometerCors import *
11from l2q import *
12from combineMulti import *
13from _procedures import *
14#from mantidsimple import *  # Old API
15from mantid.simpleapi import *  # New API
16from mantid.api import WorkspaceGroup
17import math
18import re
19
20
21def quick(run, theta=0, pointdet=1,roi=[0,0], db=[0,0], trans='', polcorr=0, usemon=-1,outputType='pd', debug=0):
22        '''
23        call signature(s)::
24
25        x=quick(RunNumber)
26        x=quick(RunNumber, roi=[0,0], db=[0,0], trans=0, outputType='pd')
27        x=quick(RunNumber,[1,10])
28        x=quick(RunNumber,[1,10],[20,40])
29        x=quick(RunNumber, trans=2568)
30        x=quick(RunNumber, trans='SomeSavedWorkspaceName')
31               
32        Reduces a ISIS  raw or nexus file created on one of the reflectometers applying
33        only a minimum ammount of corrections. The data is left in terms of lambda.
34       
35        Required arguments
36        =========   =====================================================================
37        RunNumber       Either an ISIS run number when the paths are set up correctly or
38                        the full path and filename if an ISIS raw file.
39        =========   =====================================================================
40
41        Optional keyword arguments:     
42        =========   =====================================================================
43        Keyword         Description
44        =========   =====================================================================
45        roi             Region of interest marking the extent of the reflected beam.
46                        default [0,0]
47        db              Region of interest marking the extent of the direct beam.
48                        default [0,0]
49        trans           transmission run number or saved workspace. The default is 0 (No
50                        transmission run).  trans=-1 will supress the division of the
51                        detector by the monitor.
52        polcorr         polarisation correction, 0=no correction (unpolarised run)
53        usemon          monitor to be used for normalisation (-1 is default from IDF)
54        outputType      'pd' = point detector (Default), 'md'=  Multidetector   Will use
55                        this to build the equivalent of gd in the old matlab code but
56                        keep all of the simple detector processing in one well organized
57                        function.  This should not be used by the average user.
58        =========   =====================================================================
59
60        Outputs:
61        =========   =====================================================================
62        x               Either a single mantid workspace or worspace group or an array
63                        of them.
64        =========   =====================================================================
65
66        Working examples:
67        >>> # reduce a data set with the default parameters
68        >>> x=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003014.raw")
69
70        >>> # reduce a data set with a transmission run
71        >>> t=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003010.raw")
72        >>> x=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003014.raw", trans=t)
73       
74        >>> # reduce a data set using the multidetector and output a single reflectivity
75        >>> # where the reflected beam is between channel 121 and 130.
76        >>> x=quick(/Users/trc/Dropbox/Work/PolrefTest/POLREF00003014.raw", [121,130])
77
78
79        Also see: pol
80
81        ToDo:
82                1) code for the transmisson DONE!
83                2) Similar to the genie on polref add extraction from the multidetector
84                3) need to make the variables stored in the frame work contain the run number. DONE!
85               
86        '''
87
88        ''' Notes for developers:
89
90                Naming conventions for workspaces which live in the mantid framework are as follows:
91
92                        It's nearly random.  this needs to be fixed so that name clashes do not occur. 
93                        May try adding a pair of underscores to the front of the name.
94
95        '''
96       
97       
98        [I0MonitorIndex, MultiDetectorStart, nHist] = toLam(run,'',pointdet=1) #creates wTof = "_W" + name
99        inst = groupGet("_W",'inst')
100        # Some beamline constants from IDF
101        intmin = inst.getNumberParameter('MonitorIntegralMin')[0]
102        intmax = inst.getNumberParameter('MonitorIntegralMax')[0]
103        print I0MonitorIndex
104        print nHist
105        if (nHist > 5 and not(pointdet)):
106                # Proccess Multi-Detector; assume MD goes to the end:
107                # if roi or db are given in the function then sum over the apropriate channels
108                print "This is a multidetector run."
109                try:
110                        CropWorkspace(InputWorkspace="_D",OutputWorkspace="_DM",StartWorkspaceIndex=MultiDetectorStart)
111                        RebinToWorkspace(WorkspaceToRebin="_M",WorkspaceToMatch="_DM",OutputWorkspace="_M_M")
112                        CropWorkspace(InputWorkspace="_M_M",OutputWorkspace="_I0M",StartWorkspaceIndex=I0MonitorIndex)
113                        RebinToWorkspace(WorkspaceToRebin="_I0M",WorkspaceToMatch="_DM",OutputWorkspace="_I0M")
114                        Divide(LHSWorkspace="_DM",RHSWorkspace="_I0M",OutputWorkspace="IvsLam")
115                        if (roi != [0,0]) :
116                                SumSpectra(InputWorkspace="IvsLam",OutputWorkspace="DMR",StartWorkspaceIndex=roi[0], EndWorkspaceIndex=roi[1])
117                                ReflectedBeam=mtd['DMR']
118                        if (db != [0,0]) :
119                                SumSpectra(InputWorkspace="_DM",OutputWorkspace="_DMD",StartWorkspaceIndex=db[0], EndWorkspaceIndex=db[1])
120                                DirectBeam=mtd['_DMD']
121                except SystemExit:
122                        print "Point-Detector only run."
123                RunNumber = groupGet('IvsLam','samp','run_number')
124                if (theta):
125                        IvsQ = l2q(mtd['DMR'], 'linear-detector', theta)
126                else:
127                        ConvertUnits(InputWorkspace='DMR',OutputWorkspace="IvsQ",Target="MomentumTransfer")
128                               
129        # Single Detector processing-------------------------------------------------------------
130        else:
131                print "This is a Point-Detector run."
132                # handle transmission runs
133                # process the point detector reflectivity 
134                RebinToWorkspace(WorkspaceToRebin="_M",WorkspaceToMatch="_DP",OutputWorkspace="_M_P")
135                CropWorkspace(InputWorkspace="_M_P",OutputWorkspace="_I0P",StartWorkspaceIndex=I0MonitorIndex,EndWorkspaceIndex=I0MonitorIndex)
136                Scale(InputWorkspace="_DP",OutputWorkspace="IvsLam",Factor=1)
137                #Divide(LHSWorkspace="_DP",RHSWorkspace="_I0P",OutputWorkspace="IvsLam")
138                #  Normalise by good frames
139                GoodFrames = groupGet('IvsLam','samp','goodfrm')
140                print "run frames: ", GoodFrames
141                if (run=='0'):
142                        RunNumber = '0'
143                else:
144                        RunNumber = groupGet('IvsLam','samp','run_number') 
145                #mantid['IvsLam'].getRun().getLogData("goodfrm").value
146                #Scale('IvsLam','IvsLam',GoodFrames**-1,'Multiply')
147                #IvsLam = mantid['IvsLam']*GoodFrames**-1
148                if (trans==''):
149                        #monitor2Eff('M')  # This doesn't seem to work.
150                        #heliumDetectorEff('DP')  # point detector  #Nor does this.
151                        # Multidetector   (Flood)   TODO               
152                        print "No transmission file. Trying default exponential/polynomial correction..."
153                        inst=groupGet('_DP','inst')
154                        corrType=inst.getStringParameter('correction')[0]
155                        if (corrType=='polynomial'):
156                                pString=inst.getStringParameter('polystring')
157                                print pString
158                                if len(pString):
159                                        PolynomialCorrection(InputWorkspace='_DP',OutputWorkspace='IvsLam',Coefficients=pString[0],Operation='Divide')
160                                else:
161                                        print "No polynomial coefficients in IDF. Using monitor spectrum with no corrections."
162                        elif (corrType=='exponential'):
163                                c0=inst.getNumberParameter('C0')
164                                c1=inst.getNumberParameter('C1')
165                                print "Exponential parameters: ", c0[0], c1[0]
166                                if len(c0):
167                                        ExponentialCorrection(InputWorkspace='_DP',OutputWorkspace='IvsLam',C0=c0[0],C1=c1[0],Operation='Divide')
168                                # normalise by monitor spectrum
169                                # RebinToWorkspace(WorkspaceToRebin="_M",WorkspaceToMatch="_DP",OutputWorkspace="_M_M")
170                                # CropWorkspace(InputWorkspace="M_M",OutputWorkspace="I0M",StartWorkspaceIndex=I0MonitorIndex)
171                                # Divide(LHSWorkspace="DM",RHSWorkspace="I0M",OutputWorkspace="RM")
172                        Divide(LHSWorkspace="IvsLam",RHSWorkspace="_I0P",OutputWorkspace="IvsLam")
173                else: # we have a transmission run
174                        Integration(InputWorkspace="_I0P",OutputWorkspace="_monInt",RangeLower=str(intmin),RangeUpper=str(intmax))
175                        #scaling=1/mantid.getMatrixWorkspace('_monInt').dataY(0)[0]
176                        Divide(LHSWorkspace="_DP",RHSWorkspace="_monInt",OutputWorkspace="IvsLam")
177                        ##Divide(LHSWorkspace="_DP",RHSWorkspace="_I0P",OutputWorkspace="IvsLam")
178                        names = mtd.getObjectNames()
179                        if trans in names:
180                                ##Divide(LHSWorkspace="_DP",RHSWorkspace="_I0P",OutputWorkspace="IvsLam")
181                                RebinToWorkspace(WorkspaceToRebin=trans,WorkspaceToMatch="IvsLam",OutputWorkspace=trans)
182                                ##IvsLam = mantid['IvsLam']*GoodFrames**-1
183                                Divide(LHSWorkspace="IvsLam",RHSWorkspace=trans,OutputWorkspace="IvsLam")
184                        else:
185                                transCorr(trans)
186                # Need to process the optional args to see what needs to be output and what division needs to be made
187                if (polcorr==1):
188                        doPNR('IvsLam')
189                elif (polcorr==2):
190                        doPA('IvsLam')
191                elif (polcorr==3):
192                        print "Sorry - no other correction is defined yet!"
193                # Convert to I vs Q
194                # check if detector in direct beam
195                if (theta == 0 or theta == ''):
196                        if (theta == ''):
197                                theta = 0
198                        print "given theta = ",theta
199                        inst = groupGet('IvsLam','inst')
200                        detLocation=inst.getComponentByName('point-detector').getPos()
201                        sampleLocation=inst.getComponentByName('some-surface-holder').getPos()
202                        detLocation=inst.getComponentByName('point-detector').getPos()
203                        sample2detector=detLocation-sampleLocation    # metres
204                        source=inst.getSource()
205                        beamPos = sampleLocation - source.getPos()
206                        PI = 3.1415926535
207                        theta = inst.getComponentByName('point-detector').getTwoTheta(sampleLocation, beamPos)*180.0/PI/2.0
208                        print "Det location: ", detLocation, "Calculated theta = ",theta
209                        if detLocation.getY() == 0:  # detector is not in correct place
210                                print "det at 0"
211                                # Load corresponding NeXuS file
212                                runno = '_' + str(run)
213                                #templist.append(runno)
214                                if type(run)==type(int()):
215                                        LoadNexus(Filename=run,OutputWorkspace=runno)
216                                else:
217                                        LoadNexus(Filename=run.replace("raw","nxs",1),OutputWorkspace=runno)
218                                # Get detector angle theta from NeXuS
219                                theta = groupGet(runno,'samp','theta')
220                                print 'Nexus file theta =', theta
221                                IvsQ = l2q(mtd['IvsLam'], 'point-detector', theta)
222                        else:
223                                ConvertUnits(InputWorkspace='IvsLam',OutputWorkspace="IvsQ",Target="MomentumTransfer")
224                       
225                else:
226                        theta = float(theta)
227                        IvsQ = l2q(mtd['IvsLam'], 'point-detector', theta)             
228       
229        RenameWorkspace(InputWorkspace='IvsLam',OutputWorkspace=RunNumber+'_IvsLam')
230        RenameWorkspace(InputWorkspace='IvsQ',OutputWorkspace=RunNumber+'_IvsQ')
231               
232        # delete all temporary workspaces unless in debug mode (debug=1)
233   
234        if debug == 0:
235                cleanup()
236        return  mtd[RunNumber+'_IvsLam'], mtd[RunNumber+'_IvsQ'], theta
237
238
239def doPNR(wksp):
240        inst = groupGet("_W",'inst')
241        # Some beamline constants from IDF
242        crho = inst.getStringParameter('crho')
243        calpha = inst.getStringParameter('calpha')
244        cAp = inst.getStringParameter('cAp')
245        cPp = inst.getStringParameter('cPp')
246        nrPNRCorrection(wksp,crho[0],calpha[0],cAp[0],cPp[0])
247       
248
249def doPA(wksp):
250        inst = groupGet("_W",'inst')
251        # Some beamline constants from IDF
252        crho = inst.getStringParameter('crho')
253        calpha = inst.getStringParameter('calpha')
254        cAp = inst.getStringParameter('cAp')
255        cPp = inst.getStringParameter('cPp')
256        nrPACorrection(wksp,crho[0],calpha[0],cAp[0],cPp[0])
257
258
259def nrPNRCorrection(Wksp,crho,calpha,cAp,cPp):
260
261        # Constants Based on Runs 18350+18355 and 18351+18356 analyser theta at -0.1deg
262        # 2 RF Flippers as the polarising system
263        # crho=[1.006831,-0.011467,0.002244,-0.000095]
264        # calpha=[1.017526,-0.017183,0.003136,-0.000140]
265        # cAp=[0.917940,0.038265,-0.006645,0.000282]
266        # cPp=[0.972762,0.001828,-0.000261,0.0]
267        print "Performing PNR correction with parameters from IDF..."
268        CloneWorkspace(Wksp,OutputWorkspace='_'+Wksp+'_uncorrected')
269        if (mantid[Wksp].size()!=2):
270                print "PNR correction works only with exactly 2 periods!"
271        else:
272                Ip = mtd[Wksp+'_1']
273                Ia = mtd[Wksp+'_2']
274
275                CloneWorkspace(Ip,OutputWorkspace="PCalpha")
276                CropWorkspace(InputWorkspace="PCalpha",OutputWorkspace="PCalpha",StartWorkspaceIndex="0",EndWorkspaceIndex="0")
277                PCalpha=(mtd['PCalpha']*0.0)+1.0
278                alpha=mtd['PCalpha']
279                # a1=alpha.readY(0)
280                # for i in range(0,len(a1)):
281                        # alpha.dataY(0)[i]=0.0
282                        # alpha.dataE(0)[i]=0.0
283                CloneWorkspace("PCalpha",OutputWorkspace="PCrho")
284                CloneWorkspace("PCalpha",OutputWorkspace="PCAp")
285                CloneWorkspace("PCalpha",OutputWorkspace="PCPp")
286                rho=mtd['PCrho']
287                Ap=mtd['PCAp']
288                Pp=mtd['PCPp']
289                # for i in range(0,len(a1)):
290                        # x=(alpha.dataX(0)[i]+alpha.dataX(0)[i])/2.0
291                        # for j in range(0,4):
292                                # alpha.dataY(0)[i]=alpha.dataY(0)[i]+calpha[j]*x**j
293                                # rho.dataY(0)[i]=rho.dataY(0)[i]+crho[j]*x**j
294                                # Ap.dataY(0)[i]=Ap.dataY(0)[i]+cAp[j]*x**j
295                                # Pp.dataY(0)[i]=Pp.dataY(0)[i]+cPp[j]*x**j
296                PolynomialCorrection(InputWorkspace="PCalpha",OutputWorkspace="PCalpha",Coefficients=calpha,Operation="Multiply")
297                PolynomialCorrection(InputWorkspace="PCrho",OutputWorkspace="PCrho",Coefficients=crho,Operation="Multiply")
298                PolynomialCorrection(InputWorkspace="PCAp",OutputWorkspace="PCAp",Coefficients=cAp,Operation="Multiply")
299                PolynomialCorrection(InputWorkspace="PCPp",OutputWorkspace="PCPp",Coefficients=cPp,Operation="Multiply")
300                D=Pp*(1.0+rho)
301                nIp=(Ip*(rho*Pp+1.0)+Ia*(Pp-1.0))/D
302                nIa=(Ip*(rho*Pp-1.0)+Ia*(Pp+1.0))/D
303                RenameWorkspace(nIp,OutputWorkspace=str(Ip)+"corr")
304                RenameWorkspace(nIa,OutputWorkspace=str(Ia)+"corr")
305               
306                GroupWorkspaces(str(Ip)+"corr"+','+str(Ia)+"corr",OutputWorkspace=Wksp)
307                #CloneWorkspace=(InputWorkspace="gr",OutputWorkspace=Wksp)
308                iwksp=mantid.getWorkspaceNames()
309                list=[str(Ip),str(Ia),"PCalpha","PCrho","PCAp","PCPp","1_p"]
310                for i in range(len(iwksp)):
311                        for j in list:
312                                lname=len(j)
313                                if iwksp[i] [0:lname+1] == j+"_":
314                                        mantid.deleteWorkspace(iwksp[i])
315                mantid.deleteWorkspace("PCalpha")
316                mantid.deleteWorkspace("PCrho")
317                mantid.deleteWorkspace("PCAp")
318                mantid.deleteWorkspace("PCPp")
319                mantid.deleteWorkspace("D")
320                mantid.deleteWorkspace(str(Ip)+"corr")
321                mantid.deleteWorkspace(str(Ia)+"corr")
322#
323#       
324def nrPACorrection(Wksp,crho,calpha,cAp,cPp):#UpUpWksp,UpDownWksp,DownUpWksp,DownDownWksp):
325#       crho=[0.941893,0.0234006,-0.00210536,0.0]
326#       calpha=[0.945088,0.0242861,-0.00213624,0.0]
327#       cAp=[1.00079,-0.0186778,0.00131546,0.0]
328#       cPp=[1.01649,-0.0228172,0.00214626,0.0]
329        # Constants Based on Runs 18350+18355 and 18351+18356 analyser theta at -0.1deg
330        # 2 RF Flippers as the polarising system
331        # Ipa and Iap appear to be swapped in the sequence on CRISP 4 perido data!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
332        print "Performing PA correction with parameters from IDF..."
333        CloneWorkspace(Wksp,OutputWorkspace='_'+Wksp+'_uncorrected')
334        if (mantid[Wksp].size()!=4):
335                print "PNR correction works only with exactly 4 periods (uu,ud,du,dd)!"
336        else:
337                Ipp = mtd[Wksp+'_1']
338                Ipa = mtd[Wksp+'_2']
339                Iap = mtd[Wksp+'_3']
340                Iaa = mtd[Wksp+'_4']
341               
342                CloneWorkspace(Ipp,OutputWorkspace="PCalpha")
343                CropWorkspace(InputWorkspace="PCalpha",OutputWorkspace="PCalpha",StartWorkspaceIndex="0",EndWorkspaceIndex="0")
344                PCalpha=(mtd['PCalpha']*0.0)+1.0
345                alpha=mtd['PCalpha']
346                # a1=alpha.readY(0)
347                # for i in range(0,len(a1)):
348                        # alpha.dataY(0)[i]=0.0
349                        # alpha.dataE(0)[i]=0.0
350                CloneWorkspace("PCalpha",OutputWorkspace="PCrho")
351                CloneWorkspace("PCalpha",OutputWorkspace="PCAp")
352                CloneWorkspace("PCalpha",OutputWorkspace="PCPp")
353                rho=mtd['PCrho']
354                Ap=mtd['PCAp']
355                Pp=mtd['PCPp']
356                # for i in range(0,len(a1)):
357                        # x=(alpha.dataX(0)[i]+alpha.dataX(0)[i])/2.0
358                        # for j in range(0,4):
359                                # alpha.dataY(0)[i]=alpha.dataY(0)[i]+calpha[j]*x**j
360                                # rho.dataY(0)[i]=rho.dataY(0)[i]+crho[j]*x**j
361                                # Ap.dataY(0)[i]=Ap.dataY(0)[i]+cAp[j]*x**j
362                                # Pp.dataY(0)[i]=Pp.dataY(0)[i]+cPp[j]*x**j
363                # Use the polynomial corretion fn instead
364                PolynomialCorrection(InputWorkspace="PCalpha",OutputWorkspace="PCalpha",Coefficients=calpha,Operation="Multiply")
365                PolynomialCorrection(InputWorkspace="PCrho",OutputWorkspace="PCrho",Coefficients=crho,Operation="Multiply")
366                PolynomialCorrection(InputWorkspace="PCAp",OutputWorkspace="PCAp",Coefficients=cAp,Operation="Multiply")
367                PolynomialCorrection(InputWorkspace="PCPp",OutputWorkspace="PCPp",Coefficients=cPp,Operation="Multiply")
368               
369                A0 = (Iaa * Pp * Ap) + (Ap * Ipa * rho * Pp) + (Ap * Iap * Pp * alpha) + (Ipp * Ap * alpha * rho * Pp)
370                A1 = Pp * Iaa
371                A2 = Pp * Iap
372                A3 = Ap * Iaa
373                A4 = Ap * Ipa
374                A5 = Ap * alpha * Ipp
375                A6 = Ap * alpha * Iap
376                A7 = Pp * rho  * Ipp
377                A8 = Pp * rho  * Ipa
378                D = Pp * Ap *( 1.0 + rho + alpha + (rho * alpha) ) 
379                nIpp = (A0 - A1 + A2 - A3 + A4 + A5 - A6 + A7 - A8 + Ipp + Iaa - Ipa - Iap) / D
380                nIaa = (A0 + A1 - A2 + A3 - A4 - A5 + A6 - A7 + A8 + Ipp + Iaa - Ipa - Iap) / D
381                nIpa = (A0 - A1 + A2 + A3 - A4 - A5 + A6 + A7 - A8 - Ipp - Iaa + Ipa + Iap) / D
382                nIap = (A0 + A1 - A2 - A3 + A4 + A5 - A6 - A7 + A8 - Ipp - Iaa + Ipa + Iap) / D
383                RenameWorkspace(nIpp,OutputWorkspace=str(Ipp)+"corr")
384                RenameWorkspace(nIpa,OutputWorkspace=str(Ipa)+"corr")
385                RenameWorkspace(nIap,OutputWorkspace=str(Iap)+"corr")
386                RenameWorkspace(nIaa,OutputWorkspace=str(Iaa)+"corr")
387                ReplaceSpecialValues(str(Ipp)+"corr",OutputWorkspace=str(Ipp)+"corr",NaNValue="0.0",NaNError="0.0",InfinityValue="0.0",InfinityError="0.0")
388                ReplaceSpecialValues(str(Ipp)+"corr",OutputWorkspace=str(Ipp)+"corr",NaNValue="0.0",NaNError="0.0",InfinityValue="0.0",InfinityError="0.0")
389                ReplaceSpecialValues(str(Ipp)+"corr",OutputWorkspace=str(Ipp)+"corr",NaNValue="0.0",NaNError="0.0",InfinityValue="0.0",InfinityError="0.0")
390                ReplaceSpecialValues(str(Ipp)+"corr",OutputWorkspace=str(Ipp)+"corr",NaNValue="0.0",NaNError="0.0",InfinityValue="0.0",InfinityError="0.0")
391                iwksp=mantid.getWorkspaceNames()
392                list=[str(Ipp),str(Ipa),str(Iap),str(Iaa),"PCalpha","PCrho","PCAp","PCPp","1_p"]
393                for i in range(len(iwksp)):
394                        for j in list:
395                                lname=len(j)
396                                if iwksp[i] [0:lname+1] == j+"_":
397                                        mantid.deleteWorkspace(iwksp[i])
398                mantid.deleteWorkspace("PCalpha")
399                mantid.deleteWorkspace("PCrho")
400                mantid.deleteWorkspace("PCAp")
401                mantid.deleteWorkspace("PCPp")
402                mantid.deleteWorkspace('A0')
403                mantid.deleteWorkspace('A1')
404                mantid.deleteWorkspace('A2')
405                mantid.deleteWorkspace('A3')
406                mantid.deleteWorkspace('A4')
407                mantid.deleteWorkspace('A5')
408                mantid.deleteWorkspace('A6')
409                mantid.deleteWorkspace('A7')
410                mantid.deleteWorkspace('A8')
411                mantid.deleteWorkspace('D')
412#
413
414def transCorr(transrun):
415        inst = groupGet("_W",'inst')
416        # Some beamline constants from IDF
417        intmin = inst.getNumberParameter('MonitorIntegralMin')[0]
418        intmax = inst.getNumberParameter('MonitorIntegralMax')[0]
419        if ',' in transrun:
420                slam = transrun.split(',')[0]
421                llam = transrun.split(',')[1]
422                print "Transmission runs: ", transrun
423                [I0MonitorIndex, MultiDetectorStart, nHist] = toLam(slam,'_'+slam)
424                CropWorkspace(InputWorkspace="_D_"+slam,OutputWorkspace="_D_"+slam,StartWorkspaceIndex=0,EndWorkspaceIndex=0)
425                [I0MonitorIndex, MultiDetectorStart, nHist] = toLam(llam,'_'+llam)
426                CropWorkspace(InputWorkspace="_D_"+llam,OutputWorkspace="_D_"+llam,StartWorkspaceIndex=0,EndWorkspaceIndex=0)
427               
428                RebinToWorkspace(WorkspaceToRebin="_M_"+llam,WorkspaceToMatch="_DP_"+llam,OutputWorkspace="_M_P_"+llam)
429                CropWorkspace(InputWorkspace="_M_P_"+llam,OutputWorkspace="_I0P_"+llam,StartWorkspaceIndex=I0MonitorIndex)
430               
431                #Normalise by monitor integral
432                inst = groupGet('_D_'+slam,'inst')
433                # Some beamline constants from IDF
434                intmin = inst.getNumberParameter('MonitorIntegralMin')[0]
435                intmax = inst.getNumberParameter('MonitorIntegralMax')[0]
436                Integration(InputWorkspace="_I0P_"+llam,OutputWorkspace="_monInt_TRANS",RangeLower=str(intmin),RangeUpper=str(intmax))
437                Divide(LHSWorkspace="_DP_"+llam,RHSWorkspace="_monInt_TRANS",OutputWorkspace="_D_"+llam)
438                #scaling=1/mantid.getMatrixWorkspace('_monInt_TRANS').dataY(0)[0]
439                #Scale(InputWorkspace="_DP_"+llam,OutputWorkspace="_D_"+llam,Factor=scaling,Operation="Multiply")
440               
441                # same for short wavelength run slam:
442                RebinToWorkspace(WorkspaceToRebin="_M_"+slam,WorkspaceToMatch="_DP_"+slam,OutputWorkspace="_M_P_"+slam)
443                CropWorkspace(InputWorkspace="_M_P_"+slam,OutputWorkspace="_I0P_"+slam,StartWorkspaceIndex=I0MonitorIndex)
444
445                #Normalise by monitor integral
446                inst = groupGet('_D_'+llam,'inst')
447                # Some beamline constants from IDF
448                intmin = inst.getNumberParameter('MonitorIntegralMin')[0]
449                intmax = inst.getNumberParameter('MonitorIntegralMax')[0]
450                Integration(InputWorkspace="_I0P_"+slam,OutputWorkspace="_monInt_TRANS",RangeLower=str(intmin),RangeUpper=str(intmax))
451                #scaling=1/mantid.getMatrixWorkspace('_monInt_TRANS').dataY(0)[0]
452                Divide(LHSWorkspace="_DP_"+slam,RHSWorkspace="_monInt_TRANS",OutputWorkspace="_D_"+slam)
453                #Scale(InputWorkspace="_DP_"+slam,OutputWorkspace="_D_"+slam,Factor=scaling,Operation="Multiply")
454               
455                #Divide(LHSWorkspace="_DP_"+slam,RHSWorkspace="_I0P_"+slam,OutputWorkspace="_D_"+slam)
456
457                [transr, sf] = combine2("_D_"+slam,"_D_"+llam,"_DP_TRANS",10.0,12.0,1.5,17.0,0.02,scalehigh=1)
458                #[wlam, wq, th] = quick(runno,angle,trans='_transcomb')
459        else:
460                [I0MonitorIndex, MultiDetectorStart, nHist] = toLam(transrun,'_TRANS')
461                RebinToWorkspace(WorkspaceToRebin="_M_TRANS",WorkspaceToMatch="_DP_TRANS",OutputWorkspace="_M_P_TRANS")
462                CropWorkspace(InputWorkspace="_M_P_TRANS",OutputWorkspace="_I0P_TRANS",StartWorkspaceIndex=I0MonitorIndex)
463
464                #Normalise by monitor integral
465                Integration(InputWorkspace="_I0P_TRANS",OutputWorkspace="_monInt_TRANS",RangeLower=str(intmin),RangeUpper=str(intmax))
466                Divide(LHSWorkspace="_DP_TRANS",RHSWorkspace="_monInt_TRANS",OutputWorkspace="_DP_TRANS")
467                #scaling=1/mantid.getMatrixWorkspace('_monInt_TRANS').dataY(0)[0]
468                #print "SCALING:",scaling
469                #Scale(InputWorkspace="_I0P_TRANS",OutputWorkspace=str(transrun)+"_IvsLam_TRANS",Factor=scaling,Operation="Multiply")
470                #Scale(InputWorkspace="_DP_TRANS",OutputWorkspace="_DP_TRANS",Factor=scaling,Operation="Multiply")
471               
472        #got sometimes very slight binning diferences, so do this again:
473        RebinToWorkspace(WorkspaceToRebin='_DP_TRANS',WorkspaceToMatch="IvsLam",OutputWorkspace=str(transrun)+'_IvsLam_TRANS')
474        if isinstance(mtd["_DP_TRANS"], WorkspaceGroup):
475                Divide(LHSWorkspace="IvsLam",RHSWorkspace=str(transrun)+"_IvsLam_TRANS_1",OutputWorkspace="IvsLam")
476        else:
477                Divide(LHSWorkspace="IvsLam",RHSWorkspace=str(transrun)+"_IvsLam_TRANS",OutputWorkspace="IvsLam")
478
479
480def cleanup():
481        names = mtd.getObjectNames()
482        for name in names:
483                if re.search("^_", name):
484                        DeleteWorkspace(name)
485               
486def coAdd(run,name):
487        names = mtd.getObjectNames()
488        wTof = "_W" + name              # main workspace in time-of-flight     
489        if run in names:
490                RenameWorkspace(InputWorkspace=run,OutputWorkspace=wTof)
491        else:
492
493                currentInstrument=config['default.instrument']
494                runlist = []
495                l1 = run.split(',')
496                for subs in l1:
497                        l2 = subs.split(':')
498                        for l3 in l2:
499                                runlist.append(l3)
500                print "Adding: ", runlist
501                currentInstrument=currentInstrument.upper()
502                if (runlist[0]=='0'): #DAE/current run
503                        StartLiveData(Instrument=currentInstrument,UpdateEvery='0',Outputworkspace='_sum')
504                        #LoadLiveData(currentInstrument,OutputWorkspace='_sum')
505                        #LoadDAE(DAEname='ndx'+mantid.settings['default.instrument'],OutputWorkspace='_sum')
506                else:
507                        Load(Filename=runlist[0],OutputWorkspace='_sum')#,LoadLogFiles="1")
508                       
509                for i in range(len(runlist)-1):
510                        if (runlist[i+1]=='0'): #DAE/current run
511                                StartLiveData(Instrument=currentInstrument,UpdateEvery='0',Outputworkspace='_w2')
512                                #LoadLiveData(currentInstrument,OutputWorkspace='_w2')
513                                #LoadDAE(DAEname='ndx'+mantid.settings['default.instrument'],OutputWorkspace='_w2')
514                        else:
515                                Load(Filename=runlist[i+1],OutputWorkspace='_w2')#,LoadLogFiles="1")
516                               
517                        Plus(LHSWorkspace='_sum',RHSWorkspace='_w2',OutputWorkspace='_sum')
518
519                RenameWorkspace(InputWorkspace='_sum',OutputWorkspace=wTof)
520
521       
522
523def toLam(run, name, pointdet=1):
524        '''
525        toLam splits a given run into monitor and detector spectra and
526        converts these to wavelength
527        '''
528        # some naming for convenience
529        wTof = "_W" + name              # main workspace in time-of-flight
530        monInLam = "_M" + name          # monitor spectra vs. wavelength
531        detInLam = "_D" + name          # detector spectra vs. wavelength
532        pDet = "_DP" + name                     # point-detector only vs. wavelength
533        mDet = "_DM" + name                     # multi-detector only vs. wavelength
534       
535
536        # add multiple workspaces, if given
537        coAdd(run,name)
538
539        inst = groupGet(wTof,'inst')
540        # Some beamline constants from IDF
541       
542        bgmin = inst.getNumberParameter('MonitorBackgroundMin')[0]
543        bgmax = inst.getNumberParameter('MonitorBackgroundMax')[0]
544        MonitorBackground = [bgmin,bgmax]
545        intmin = inst.getNumberParameter('MonitorIntegralMin')[0]
546        intmax = inst.getNumberParameter('MonitorIntegralMax')[0]
547        MonitorsToCorrect = [int(inst.getNumberParameter('MonitorsToCorrect')[0])]
548        # Note: Since we are removing the monitors in the load raw command they are not counted here.
549        PointDetectorStart = int(inst.getNumberParameter('PointDetectorStart')[0])
550        PointDetectorStop = int(inst.getNumberParameter('PointDetectorStop')[0])
551        MultiDetectorStart = int(inst.getNumberParameter('MultiDetectorStart')[0])
552        I0MonitorIndex = int(inst.getNumberParameter('I0MonitorIndex')[0])     
553       
554        # Get usable wavelength range
555        LambdaMin = float(inst.getNumberParameter('LambdaMin')[0])
556        LambdaMax = float(inst.getNumberParameter('LambdaMax')[0])
557                       
558        # Convert spectra from TOF to wavelength
559        ConvertUnits(InputWorkspace=wTof,OutputWorkspace=wTof+"_lam",Target="Wavelength", AlignBins='1')
560        # Separate detector an monitor spectra manually
561        CropWorkspace(InputWorkspace=wTof+"_lam",OutputWorkspace=monInLam,StartWorkspaceIndex='0',EndWorkspaceIndex=PointDetectorStart-1)
562        CropWorkspace(InputWorkspace=wTof+"_lam",OutputWorkspace=detInLam,XMin=LambdaMin,XMax=LambdaMax,StartWorkspaceIndex=PointDetectorStart)
563        # Subtract flat background from fit in range given from Instrument Def/Par File
564        CalculateFlatBackground(InputWorkspace=monInLam,OutputWorkspace=monInLam,WorkspaceIndexList=MonitorsToCorrect,StartX=MonitorBackground[0],EndX=MonitorBackground[1])
565       
566        # Is it a multidetector run?
567        nHist = groupGet(wTof+"_lam",'wksp','')
568        if (nHist<6 or (nHist>5 and pointdet)):
569                CropWorkspace(InputWorkspace=wTof+"_lam",OutputWorkspace=pDet,XMin=LambdaMin,XMax=LambdaMax,StartWorkspaceIndex=PointDetectorStart,EndWorkspaceIndex=PointDetectorStop)
570        else:
571                CropWorkspace(InputWorkspace=wTof+"_lam",OutputWorkspace=mDet,XMin=LambdaMin,XMax=LambdaMax,StartWorkspaceIndex=MultiDetectorStart)
572
573       
574       
575       
576        return I0MonitorIndex, MultiDetectorStart, nHist
577       
578def groupGet(wksp,whattoget,field=''):
579        '''
580        returns information about instrument or sample details for a given workspace wksp,
581        also if the workspace is a group (info from first group element)
582        '''
583        if (whattoget == 'inst'):
584                if isinstance(mtd[wksp], WorkspaceGroup):
585                        return mtd[wksp+'_1'].getInstrument()
586                else:
587                        return mtd[wksp].getInstrument()
588                       
589        elif (whattoget == 'samp' and field != ''):
590                if isinstance(mtd[wksp], WorkspaceGroup):
591                        try:
592                                log = mtd[wksp + '_1'].getRun().getLogData(field).value
593                                if (type(log) is int or type(log) is str):
594                                        res=log
595                                else:
596                                        res = log[len(log)-1]
597                        except RuntimeError:
598                                res = 0
599                                print "Block "+field+" not found."                     
600                else:
601                        try:
602                                log = mtd[wksp].getRun().getLogData(field).value
603                                if (type(log) is int or type(log) is str):
604                                        res=log
605                                else:
606                                        res = log[len(log)-1]
607                        except RuntimeError:           
608                                res = 0
609                                print "Block "+field+" not found."
610                return res
611        elif (whattoget == 'wksp'):
612                if isinstance(mtd[wksp], WorkspaceGroup):
613                        return mtd[wksp+'_1'].getNumberHistograms()
614                else:
615                        return mtd[wksp].getNumberHistograms()
616
617
618""" =====================  Testing stuff =====================
619       
620        ToDo:
621                1) Make the test below more rigorous.
622                2) Each test should either print out Success or Failure
623                or describe in words what the tester should see on the screen. 
624                3) Each function should be accompanied by a test function
625                4) all test functions to be run in a give file should be called in doAllTests().
626                The reason for this is that the normal python if __name__ == '__main__': 
627                testing location is heavily used during script development and debugging.
628
629"""
630       
631def _testQuick():
632        config['default.instrument'] = "SURF"
633        [w1lam,w1q,th] = quick(94511,theta=0.25,trans='94504')
634        [w2lam,w2q,th] = quick(94512,theta=0.65,trans='94504')
635        [w3lam,w3q,th] = quick(94513,theta=1.5,trans='94504')
636        g1=plotSpectrum("94511_IvsQ",0)
637        g2=plotSpectrum("94512_IvsQ",0)
638        g3=plotSpectrum("94513_IvsQ",0)
639       
640        return True
641
642def  _doAllTests():
643        _testQuick()
644        return True
645
646
647if __name__ == '__main__':
648        ''' This is the debugging and testing area of the file.  The code below is run when ever the
649           script is called directly from a shell command line or the execute all menu option in mantid.
650        '''
651        #Debugging = True  # Turn the debugging on and the testing code off
652        Debugging = False # Turn the debugging off and the testing on
653       
654        if Debugging == False:
655                _doAllTests() 
656        else:    #Debugging code goes below
657                rr = quick("N:/SRF93080.raw")
658       
659
660#  x=quick(95266)