Ticket #9643: verify9643.py

File verify9643.py, 8.9 KB (added by Wenduo Zhou, 6 years ago)
Line 
1# This is the serious verification script for ticket 9643 for supporting VULCAN
2
3import math
4
5# Constants/Input
6eventfname = "/SNS/VULCAN/IPTS-10560/0/51624/NeXus/VULCAN_51624_event.nxs"
7vuloffsetfname = "pid_offset_vulcan_new.dat"
8effgeofname = "mod_difc_vulcan_new.dat"
9
10def 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
55def 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
118def 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
287if __name__ == "__main__":
288    numbanks = 6 # option: 1, 2, 6
289    main(numbanks)