1 | # This is the serious verification script for ticket 9643 for supporting VULCAN |
---|
2 | |
---|
3 | import math |
---|
4 | |
---|
5 | # Constants/Input |
---|
6 | eventfname = "/SNS/VULCAN/IPTS-10560/0/51624/NeXus/VULCAN_51624_event.nxs" |
---|
7 | vuloffsetfname = "pid_offset_vulcan_new.dat" |
---|
8 | effgeofname = "mod_difc_vulcan_new.dat" |
---|
9 | |
---|
10 | def ReadEffectiveGeometryFile(filename, throwemptybank): |
---|
11 | """ Read the effective geometry information file from VULCAN |
---|
12 | Example: mod_difc_vulcan_new.dat |
---|
13 | """ |
---|
14 | # Import file contents |
---|
15 | try: |
---|
16 | gfile = open(filename, "r") |
---|
17 | lines = gfile.readlines() |
---|
18 | gfile.close() |
---|
19 | except IOError, err: |
---|
20 | print "Geometry file %s cannot be opened or read. " % (filename) |
---|
21 | return (False, None) |
---|
22 | |
---|
23 | # Dictionary to return |
---|
24 | geomdict = {} |
---|
25 | |
---|
26 | # Import information |
---|
27 | for line in lines: |
---|
28 | line = line.strip() |
---|
29 | |
---|
30 | # reject empty line and comment line |
---|
31 | if len(line) == 0: |
---|
32 | continue |
---|
33 | terms = line.split() |
---|
34 | if terms[0].isdigit() is False: |
---|
35 | continue |
---|
36 | |
---|
37 | try: |
---|
38 | bankid = int(terms[0]) |
---|
39 | difc = float(terms[2]) |
---|
40 | refpid = int(terms[3]) |
---|
41 | twotheta = float(terms[4]) |
---|
42 | except Exception: |
---|
43 | print "Line '%s' is not conformed with standard format. " % (line) |
---|
44 | continue |
---|
45 | |
---|
46 | # reject non-used bank |
---|
47 | if throwemptybank is True and refpid == 0: |
---|
48 | continue |
---|
49 | |
---|
50 | geomdict[bankid] = {"DIFC": difc, "2Theta": twotheta, "RefPID": refpid} |
---|
51 | |
---|
52 | return (True, geomdict) |
---|
53 | |
---|
54 | |
---|
55 | def ManualDiffractionFocus(eventwsname, groupwsname, groupdict=None): |
---|
56 | """ Manually diffraction focuss the event workspace in unit of TOF |
---|
57 | Arguement: |
---|
58 | - groupdict :: Dictionary is for second level group such that to sum spectra for |
---|
59 | 2 banks and 1 bank option. |
---|
60 | """ |
---|
61 | eventws = mtd[eventwsname] |
---|
62 | groupws = mtd[groupwsname] |
---|
63 | |
---|
64 | # Check unit |
---|
65 | evunit = eventws.getAxis(0).getUnit().unitID() |
---|
66 | if evunit != "TOF": |
---|
67 | raise NotImplementedError("EventWorkspace must be TOF. It is %s now!" % (evunit)) |
---|
68 | |
---|
69 | # Sum spectrum |
---|
70 | groupspecdict = {} |
---|
71 | for groupid in xrange(1, 7): |
---|
72 | groupspecdict[groupid] = "" |
---|
73 | |
---|
74 | # Get list of spectrum |
---|
75 | for iws in xrange(groupws.getNumberHistograms()): |
---|
76 | gid = int(groupws.readY(iws)[0]) |
---|
77 | groupspecdict[gid] += "%d," % (iws) |
---|
78 | |
---|
79 | for groupid in xrange(1, 7): |
---|
80 | groupspecdict[groupid] = groupspecdict[groupid][0:-1] |
---|
81 | |
---|
82 | #for groupid in xrange(1, 7): |
---|
83 | # print "Group %d: %s" % (groupid, groupspecdict[groupid]) |
---|
84 | |
---|
85 | for groupid in xrange(1, 7): |
---|
86 | if len(groupspecdict[groupid]) == 0: |
---|
87 | print "No spectrum in group/bank %d. Skip!" % (groupid) |
---|
88 | continue |
---|
89 | outwsname = eventwsname+"_bank%d" % (groupid) |
---|
90 | SumSpectra(InputWorkspace=eventws, OutputWorkspace=outwsname, |
---|
91 | ListOfWorkspaceIndices=groupspecdict[groupid]) |
---|
92 | |
---|
93 | outws = mtd[outwsname] |
---|
94 | print "Workspace %s has %d events. " % (outwsname, outws.getNumberEvents()) |
---|
95 | |
---|
96 | return |
---|
97 | |
---|
98 | # The followed is useless |
---|
99 | # Second-level grouping |
---|
100 | if groupdict is not None: |
---|
101 | for subgroupname in groupdict.keys(): |
---|
102 | subgrouplist = groupdict[subgroupname] |
---|
103 | outputwsname = eventwsname+"_%s"%(subgroupname) |
---|
104 | if len(subgrouplist) == 0: |
---|
105 | raise NotImplementedError("Empty subgroup!") |
---|
106 | else: |
---|
107 | wsname0 = eventwsname+"_bank%d" % (subgrouplist[0]) |
---|
108 | CloneWorkspace(InputWorkspace = wsname0, OutputWorkspace = outputwsname) |
---|
109 | for gid in subgrouplist[1:]: |
---|
110 | wsnamex = eventwsname + "_bank%d" % (gid) |
---|
111 | Plus(LHSWorkspace = outputwsname, RHSWorkspace = wsnamex, OutputWorkspace = outputwsname) |
---|
112 | |
---|
113 | # ENDIF |
---|
114 | |
---|
115 | return |
---|
116 | |
---|
117 | |
---|
118 | def main(bankmode): |
---|
119 | """ verification main |
---|
120 | - bankmode :: integer (1, 2, 6) |
---|
121 | """ |
---|
122 | # Import VULCAN DIFC geometry file |
---|
123 | tup = ReadEffectiveGeometryFile(effgeofname, True) |
---|
124 | if tup[0] is False: |
---|
125 | print "Unable to read VUCLAN DIFC geometry file %s. Stop!" % (effgeofname) |
---|
126 | return |
---|
127 | |
---|
128 | geomdict = tup[1] |
---|
129 | print "Used banks: ", sorted(geomdict.keys()) |
---|
130 | |
---|
131 | # Load event data |
---|
132 | eventwsname = eventfname.split("/")[-1].split(".")[0] |
---|
133 | Load(Filename = eventfname, OutputWorkspace=eventwsname) |
---|
134 | eventws = mtd[eventwsname] |
---|
135 | |
---|
136 | # Calcualte L2 |
---|
137 | L1 = eventws.getInstrument().getSource().getPos().norm() |
---|
138 | for bankid in sorted(geomdict.keys()): |
---|
139 | difc = geomdict[bankid]["DIFC"] |
---|
140 | twotheta = geomdict[bankid]["2Theta"] |
---|
141 | l = difc/252.777/2./math.sin(twotheta*0.5*math.pi/180.) |
---|
142 | l2 = l - L1 |
---|
143 | geomdict[bankid]["L2"] = l2 |
---|
144 | print "Bank %d, L = %f, L2 = %f" % (bankid, l, l2) |
---|
145 | |
---|
146 | ################################################################################### |
---|
147 | # Load Vulcan Calibration Filename and Focus event in unit of TOF |
---|
148 | ################################################################################### |
---|
149 | # Convert the effective DIFCs and 2thetas to N-bank mode |
---|
150 | newgeomdict = {} |
---|
151 | if bankmode == 2: |
---|
152 | for bankid in [21, 22, 23]: |
---|
153 | newgeomdict[bankid] = geomdict[22] |
---|
154 | for bankid in [26, 27, 28]: |
---|
155 | newgeomdict[bankid] = geomdict[27] |
---|
156 | elif bankmode == 1: |
---|
157 | for bankid in geomdict.keys(): |
---|
158 | newgeomdict[bankid] = geomdict[22] |
---|
159 | elif bankmode == 6: |
---|
160 | newgeomdict = geomdict |
---|
161 | else: |
---|
162 | raise NotImplementedError("Bank mode %d is not supported. " % (bankmode)) |
---|
163 | |
---|
164 | # Write out the input string for LoadVulcanCalFile |
---|
165 | bankids = "" |
---|
166 | effdifcs = "" |
---|
167 | eff2thetas = "" |
---|
168 | for bankid in newgeomdict: |
---|
169 | difc = newgeomdict[bankid]["DIFC"] |
---|
170 | twotheta = newgeomdict[bankid]["2Theta"] |
---|
171 | bankids+= "%d," % (bankid) |
---|
172 | effdifcs +="%f," % (difc) |
---|
173 | eff2thetas += "%f," % (twotheta) |
---|
174 | |
---|
175 | # remove last ',' |
---|
176 | bankids = bankids[0:-1] |
---|
177 | effdifcs = effdifcs[0:-1] |
---|
178 | eff2thetas = eff2thetas[0:-1] |
---|
179 | |
---|
180 | print "Bank IDs = ", bankids |
---|
181 | print "DIFCs = ", effdifcs |
---|
182 | print "2Thetas = ", eff2thetas |
---|
183 | |
---|
184 | tofcal_eventwsname = eventwsname+"_TOF_Cal" |
---|
185 | CloneWorkspace(InputWorkspace=eventwsname, OutputWorkspace=tofcal_eventwsname) |
---|
186 | |
---|
187 | # LoadVulcanCalFile and work on an test EventWorkspace |
---|
188 | if bankmode == 2: |
---|
189 | groupmode = '2Banks' |
---|
190 | elif bankmode == 1: |
---|
191 | groupmode = '1Bank' |
---|
192 | elif bankmode == 6: |
---|
193 | groupmode = '6Modules' |
---|
194 | LoadVulcanCalFile( |
---|
195 | OffsetFilename=r'pid_offset_vulcan_new.dat', |
---|
196 | BadPixelFilename=r'bad_pids_vulcan_new_6867_7323.dat', |
---|
197 | Grouping = groupmode, |
---|
198 | BankIDs = bankids, |
---|
199 | EffectiveDIFCs = effdifcs, |
---|
200 | Effective2Thetas = eff2thetas, |
---|
201 | WorkspaceName='Vulcan_idl', |
---|
202 | EventWorkspace=tofcal_eventwsname) |
---|
203 | |
---|
204 | Rebin(InputWorkspace=tofcal_eventwsname, OutputWorkspace=tofcal_eventwsname, Params="-0.001", PreserveEvents=True) |
---|
205 | |
---|
206 | # Do diffraction focussing manually |
---|
207 | groupwsname = "Vulcan_idl_group" |
---|
208 | tofcaleventwsname = eventfname.split("/")[-1].split(".")[0] + "_TOF_Cal" |
---|
209 | |
---|
210 | if groupmode == 1: |
---|
211 | subgroupdict = {"Pack": [1,2,3,4,5,6]} |
---|
212 | elif groupmode == 2: |
---|
213 | subgroupdict = {"East": [2, 1, 3], "West": [5,4,6]} |
---|
214 | else: |
---|
215 | subgroupdict = None |
---|
216 | ManualDiffractionFocus(tofcaleventwsname, groupwsname, subgroupdict) |
---|
217 | |
---|
218 | |
---|
219 | ################################################################################### |
---|
220 | # Focus event data by using SNSPowderReduction's approach (for vervification) |
---|
221 | ################################################################################### |
---|
222 | if bankmode == 1: |
---|
223 | banks = [22] |
---|
224 | elif bankmode == 2: |
---|
225 | banks = [22, 27] |
---|
226 | elif bankmode == 6: |
---|
227 | banks = geomdict.keys() |
---|
228 | |
---|
229 | eff2thetas = "" |
---|
230 | effl2s = "" |
---|
231 | for bid in sorted(banks): |
---|
232 | twotheta = geomdict[bid]["2Theta"] |
---|
233 | l2 = geomdict[bid]["L2"] |
---|
234 | effl2s += "%f," % (l2) |
---|
235 | eff2thetas += "%f," %(twotheta) |
---|
236 | effl2s = effl2s[0:-1] |
---|
237 | eff2thetas = eff2thetas[0:-1] |
---|
238 | |
---|
239 | print "L2s: ", effl2s |
---|
240 | print "2Thetas: ", eff2thetas |
---|
241 | |
---|
242 | groupwsname = "Vulcan_idl_group" |
---|
243 | offsetswsname = "Vulcan_idl_offsets" |
---|
244 | print "Event workspace: ", eventwsname |
---|
245 | AlignDetectors(InputWorkspace = eventwsname, OutputWorkspace = eventwsname, |
---|
246 | OffsetsWorkspace = offsetswsname) |
---|
247 | DiffractionFocussing(InputWorkspace = eventwsname, OutputWorkspace = eventwsname, |
---|
248 | GroupingWorkspace = groupwsname, PreserveEvents=True) |
---|
249 | EditInstrumentGeometry(Workspace = eventwsname, |
---|
250 | L2 = effl2s, Polar = eff2thetas) |
---|
251 | ConvertUnits(InputWorkspace=eventwsname, OutputWorkspace=eventwsname, |
---|
252 | Target="TOF", Emode="Elastic") |
---|
253 | |
---|
254 | ################################################################################### |
---|
255 | # Compare |
---|
256 | ################################################################################### |
---|
257 | eventws = mtd[eventwsname] |
---|
258 | for iws in xrange(eventws.getNumberHistograms()): |
---|
259 | numevents = eventws.getEventList(iws).getNumberEvents() |
---|
260 | print "Spectrum %d: Number of events = %d" % (iws, numevents) |
---|
261 | |
---|
262 | binpar = "10000., -0.001, 50000." |
---|
263 | Rebin(InputWorkspace=eventws, OutputWorkspace=eventws, Params=binpar) |
---|
264 | |
---|
265 | |
---|
266 | for bankid in xrange(1, bankmode+1): |
---|
267 | tofcaleventws = mtd["VULCAN_51624_event_TOF_Cal_bank%d"%(bankid)] |
---|
268 | Rebin(InputWorkspace=tofcaleventws, OutputWorkspace=tofcaleventws, Params=binpar) |
---|
269 | |
---|
270 | Minus(LHSWorkspace=eventws, RHSWorkspace=tofcaleventws, AllowDifferentNumberSpectra=True, OutputWorkspace="diff_%d"%(bankid)) |
---|
271 | ConvertToMatrixWorkspace(InputWorkspace="diff_%d"%(bankid), OutputWorkspace="diff_%d"%(bankid)) |
---|
272 | |
---|
273 | # combine the difference workspaces |
---|
274 | targetdiffws = mtd["diff_%d"%(1)] |
---|
275 | for bankid in xrange(2, bankmode+1): |
---|
276 | sourcediffws = mtd["diff_%d"%(bankid)] |
---|
277 | iws = bankid-1 |
---|
278 | for i in xrange(len(sourcediffws.readY(iws))): |
---|
279 | targetdiffws.dataY(iws)[i] = sourcediffws.readY(iws)[i] |
---|
280 | DeleteWorkspace(Workspace=sourcediffws) |
---|
281 | |
---|
282 | RenameWorkspace(InputWorkspace=targetdiffws, OutputWorkspace="Diff_%s"%(eventwsname)) |
---|
283 | |
---|
284 | return |
---|
285 | |
---|
286 | |
---|
287 | if __name__ == "__main__": |
---|
288 | numbanks = 6 # option: 1, 2, 6 |
---|
289 | main(numbanks) |
---|