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 | |=============================================================================|=======| |
---|
7 | 1 80 <tab> |
---|
8 | ''' |
---|
9 | #these need to be moved into one NR folder or so |
---|
10 | #from ReflectometerCors import * |
---|
11 | from l2q import * |
---|
12 | from combineMulti import * |
---|
13 | from _procedures import * |
---|
14 | #from mantidsimple import * # Old API |
---|
15 | from mantid.simpleapi import * # New API |
---|
16 | from mantid.api import WorkspaceGroup |
---|
17 | import math |
---|
18 | import re |
---|
19 | |
---|
20 | |
---|
21 | def 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 | |
---|
239 | def 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 | |
---|
249 | def 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 | |
---|
259 | def 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 | # |
---|
324 | def 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 | |
---|
414 | def 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 | |
---|
480 | def cleanup(): |
---|
481 | names = mtd.getObjectNames() |
---|
482 | for name in names: |
---|
483 | if re.search("^_", name): |
---|
484 | DeleteWorkspace(name) |
---|
485 | |
---|
486 | def 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 | |
---|
523 | def 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 | |
---|
578 | def 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 | |
---|
631 | def _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 | |
---|
642 | def _doAllTests(): |
---|
643 | _testQuick() |
---|
644 | return True |
---|
645 | |
---|
646 | |
---|
647 | if __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) |
---|