1 | """ |
---|
2 | TXT/IGOR 2D Q Map file reader |
---|
3 | """ |
---|
4 | |
---|
5 | """ |
---|
6 | This software was developed by the University of Tennessee as part of the |
---|
7 | Distributed Data Analysis of Neutron Scattering Experiments (DANSE) |
---|
8 | project funded by the US National Science Foundation. |
---|
9 | |
---|
10 | If you use DANSE applications to do scientific research that leads to |
---|
11 | publication, we ask that you acknowledge the use of the software with the |
---|
12 | following sentence: |
---|
13 | |
---|
14 | "This work benefited from DANSE software developed under NSF award DMR-0520547." |
---|
15 | |
---|
16 | copyright 2008, University of Tennessee |
---|
17 | """ |
---|
18 | |
---|
19 | import os, sys |
---|
20 | import numpy |
---|
21 | import math, logging |
---|
22 | from DataLoader.data_info import Data2D, Detector |
---|
23 | |
---|
24 | # Look for unit converter |
---|
25 | has_converter = True |
---|
26 | try: |
---|
27 | from data_util.nxsunit import Converter |
---|
28 | except: |
---|
29 | has_converter = False |
---|
30 | |
---|
31 | class Reader: |
---|
32 | """ Simple data reader for Igor data files """ |
---|
33 | ## File type |
---|
34 | type_name = "IGOR/DAT 2D Q_map" |
---|
35 | ## Wildcards |
---|
36 | type = ["IGOR/DAT 2D file in Q_map (*.dat)|*.DAT"] |
---|
37 | ## Extension |
---|
38 | ext=['.DAT', '.dat'] |
---|
39 | |
---|
40 | def read(self,filename=None): |
---|
41 | """ Read file """ |
---|
42 | if not os.path.isfile(filename): |
---|
43 | raise ValueError, \ |
---|
44 | "Specified file %s is not a regular file" % filename |
---|
45 | |
---|
46 | # Read file |
---|
47 | f = open(filename,'r') |
---|
48 | buf = f.read() |
---|
49 | |
---|
50 | # Instantiate data object |
---|
51 | output = Data2D() |
---|
52 | output.filename = os.path.basename(filename) |
---|
53 | detector = Detector() |
---|
54 | if len(output.detector)>0: print str(output.detector[0]) |
---|
55 | output.detector.append(detector) |
---|
56 | |
---|
57 | # Get content |
---|
58 | dataStarted = False |
---|
59 | |
---|
60 | ## Defaults |
---|
61 | lines = buf.split('\n') |
---|
62 | itot = 0 |
---|
63 | x = [] |
---|
64 | y = [] |
---|
65 | |
---|
66 | ncounts = 0 |
---|
67 | |
---|
68 | wavelength = None |
---|
69 | distance = None |
---|
70 | transmission = None |
---|
71 | |
---|
72 | pixel_x = None |
---|
73 | pixel_y = None |
---|
74 | |
---|
75 | i_x = 0 |
---|
76 | i_y = -1 |
---|
77 | pixels = 0 |
---|
78 | |
---|
79 | isInfo = False |
---|
80 | isCenter = False |
---|
81 | |
---|
82 | data_conv_q = None |
---|
83 | data_conv_i = None |
---|
84 | |
---|
85 | # Set units: This is the unit assumed for Q and I in the data file. |
---|
86 | if has_converter == True and output.Q_unit != '1/A': |
---|
87 | data_conv_q = Converter('1/A') |
---|
88 | # Test it |
---|
89 | data_conv_q(1.0, output.Q_unit) |
---|
90 | |
---|
91 | if has_converter == True and output.I_unit != '1/cm': |
---|
92 | data_conv_i = Converter('1/cm') |
---|
93 | # Test it |
---|
94 | data_conv_i(1.0, output.I_unit) |
---|
95 | |
---|
96 | #Set the space for data |
---|
97 | data = numpy.zeros(0) |
---|
98 | qx_data = numpy.zeros(0) |
---|
99 | qy_data = numpy.zeros(0) |
---|
100 | q_data = numpy.zeros(0) |
---|
101 | dqx_data = numpy.zeros(0) |
---|
102 | dqy_data = numpy.zeros(0) |
---|
103 | mask = numpy.zeros(0,dtype=bool) |
---|
104 | |
---|
105 | #Read Header and 2D data |
---|
106 | for line in lines: |
---|
107 | ## Reading the header applies only to IGOR/NIST 2D q_map data files |
---|
108 | # Find setup info line |
---|
109 | if isInfo: |
---|
110 | isInfo = False |
---|
111 | line_toks = line.split() |
---|
112 | # Wavelength in Angstrom |
---|
113 | try: |
---|
114 | wavelength = float(line_toks[1]) |
---|
115 | # Units |
---|
116 | if has_converter==True and output.source.wavelength_unit != 'A': |
---|
117 | conv = Converter('A') |
---|
118 | wavelength = conv(wavelength, units=output.source.wavelength_unit) |
---|
119 | except: |
---|
120 | #Not required |
---|
121 | pass |
---|
122 | # Distance in mm |
---|
123 | try: |
---|
124 | distance = float(line_toks[3]) |
---|
125 | # Units |
---|
126 | if has_converter==True and detector.distance_unit != 'm': |
---|
127 | conv = Converter('m') |
---|
128 | distance = conv(distance, units=detector.distance_unit) |
---|
129 | except: |
---|
130 | #Not required |
---|
131 | pass |
---|
132 | |
---|
133 | # Distance in meters |
---|
134 | try: |
---|
135 | transmission = float(line_toks[4]) |
---|
136 | except: |
---|
137 | #Not required |
---|
138 | pass |
---|
139 | |
---|
140 | if line.count("LAMBDA")>0: |
---|
141 | isInfo = True |
---|
142 | |
---|
143 | # Find center info line |
---|
144 | if isCenter: |
---|
145 | isCenter = False |
---|
146 | line_toks = line.split() |
---|
147 | # Center in bin number |
---|
148 | center_x = float(line_toks[0]) |
---|
149 | center_y = float(line_toks[1]) |
---|
150 | |
---|
151 | if line.count("BCENT")>0: |
---|
152 | isCenter = True |
---|
153 | |
---|
154 | # Find data start |
---|
155 | if line.count("Data columns") or line.count("ASCII data")>0: |
---|
156 | dataStarted = True |
---|
157 | continue |
---|
158 | |
---|
159 | ## Read and get data. |
---|
160 | if dataStarted == True: |
---|
161 | line_toks = line.split() |
---|
162 | if len(line_toks) == 0: |
---|
163 | #empty line |
---|
164 | continue |
---|
165 | elif len(line_toks) < 3: |
---|
166 | #Q-map 2d require 3 or more columns of data |
---|
167 | raise ValueError,"Igor_Qmap_Reader: Can't read this file: Not a proper file format" |
---|
168 | |
---|
169 | # defaults |
---|
170 | dqx_value = None |
---|
171 | dqy_value = None |
---|
172 | qz_value = None |
---|
173 | mask_value = 0 |
---|
174 | |
---|
175 | ##Read and get data values |
---|
176 | qx_value = float(line_toks[0]) |
---|
177 | qy_value = float(line_toks[1]) |
---|
178 | value = float(line_toks[2]) |
---|
179 | |
---|
180 | try: |
---|
181 | #Read qz_value if exist on 4th column |
---|
182 | qz_value = float(line_toks[3]) |
---|
183 | except: |
---|
184 | # Found a non-float entry, skip it: Not required. |
---|
185 | pass |
---|
186 | try: |
---|
187 | #Read qz_value if exist on 5th column |
---|
188 | dqx_value = float(line_toks[4]) |
---|
189 | except: |
---|
190 | # Found a non-float entry, skip it: Not required. |
---|
191 | pass |
---|
192 | try: |
---|
193 | #Read qz_value if exist on 6th column |
---|
194 | dqy_value = float(line_toks[5]) |
---|
195 | except: |
---|
196 | # Found a non-float entry, skip it: Not required. |
---|
197 | pass |
---|
198 | try: |
---|
199 | #Read beam block mask if exist on 7th column |
---|
200 | mask_value = float(line_toks[6]) |
---|
201 | except: |
---|
202 | # Found a non-float entry, skip it |
---|
203 | pass |
---|
204 | |
---|
205 | # get data |
---|
206 | data = numpy.append(data, value) |
---|
207 | qx_data = numpy.append(qx_data, qx_value) |
---|
208 | qy_data = numpy.append(qy_data, qy_value) |
---|
209 | |
---|
210 | # optional data |
---|
211 | if dqx_value != None and numpy.isfinite(dqx_value): |
---|
212 | dqx_data = numpy.append(dqx_data, dqx_value) |
---|
213 | if dqy_value != None and numpy.isfinite(dqy_value): |
---|
214 | dqy_data = numpy.append(dqy_data, dqy_value) |
---|
215 | |
---|
216 | # default data |
---|
217 | if qz_value == None or not numpy.isfinite(qz_value): |
---|
218 | qz_value = 0 |
---|
219 | if not numpy.isfinite(mask_value): |
---|
220 | mask_value = 1 |
---|
221 | q_data = numpy.append(q_data,numpy.sqrt(qx_value**2+qy_value**2+qz_value**2)) |
---|
222 | # Note: For convenience, mask = False stands for masked, while mask = True for unmasked |
---|
223 | mask = numpy.append(mask,(mask_value>=1)) |
---|
224 | |
---|
225 | f.close() |
---|
226 | # If all mask elements are False, put all True |
---|
227 | if not mask.any(): mask[mask==False] = True |
---|
228 | |
---|
229 | # Store limits of the image in q space |
---|
230 | xmin = numpy.min(qx_data) |
---|
231 | xmax = numpy.max(qx_data) |
---|
232 | ymin = numpy.min(qy_data) |
---|
233 | ymax = numpy.max(qy_data) |
---|
234 | |
---|
235 | # units |
---|
236 | if has_converter == True and output.Q_unit != '1/A': |
---|
237 | xmin = data_conv_q(xmin, units=output.Q_unit) |
---|
238 | xmax = data_conv_q(xmax, units=output.Q_unit) |
---|
239 | ymin = data_conv_q(ymin, units=output.Q_unit) |
---|
240 | ymax = data_conv_q(ymax, units=output.Q_unit) |
---|
241 | |
---|
242 | ## calculate the range of the qx and qy_data |
---|
243 | x_size = math.fabs(xmax - xmin) |
---|
244 | y_size = math.fabs(ymax - ymin) |
---|
245 | |
---|
246 | # calculate the number of pixels in the each axes |
---|
247 | npix_y = math.floor(math.sqrt(len(data))) |
---|
248 | npix_x = math.floor(len(data)/npix_y) |
---|
249 | |
---|
250 | # calculate the size of bins |
---|
251 | xstep = x_size/(npix_x-1) |
---|
252 | ystep = y_size/(npix_y-1) |
---|
253 | |
---|
254 | # store x and y axis bin centers in q space |
---|
255 | x_bins = numpy.arange(xmin,xmax+xstep,xstep) |
---|
256 | y_bins = numpy.arange(ymin,ymax+ystep,ystep) |
---|
257 | |
---|
258 | # get the limits of q values |
---|
259 | xmin = xmin - xstep/2 |
---|
260 | xmax = xmax + xstep/2 |
---|
261 | ymin = ymin - ystep/2 |
---|
262 | ymax = ymax + ystep/2 |
---|
263 | |
---|
264 | #Store data in outputs |
---|
265 | #TODO: Check the lengths |
---|
266 | output.data = data |
---|
267 | output.err_data = numpy.sqrt(numpy.abs(data)) |
---|
268 | output.qx_data = qx_data |
---|
269 | output.qy_data = qy_data |
---|
270 | output.q_data = q_data |
---|
271 | output.mask = mask |
---|
272 | |
---|
273 | output.x_bins = x_bins |
---|
274 | output.y_bins = y_bins |
---|
275 | |
---|
276 | output.xmin = xmin |
---|
277 | output.xmax = xmax |
---|
278 | output.ymin = ymin |
---|
279 | output.ymax = ymax |
---|
280 | |
---|
281 | output.source.wavelength = wavelength |
---|
282 | |
---|
283 | # Store pixel size in mm |
---|
284 | detector.pixel_size.x = pixel_x |
---|
285 | detector.pixel_size.y = pixel_y |
---|
286 | |
---|
287 | # Store the sample to detector distance |
---|
288 | detector.distance = distance |
---|
289 | |
---|
290 | # optional data: if any of dq data == 0, do not pass to output |
---|
291 | if len(dqx_data) == len(qx_data): |
---|
292 | output.dqx_data = dqx_data |
---|
293 | if len(dqy_data) == len(qy_data): |
---|
294 | output.dqy_data = dqy_data |
---|
295 | |
---|
296 | # Units of axes |
---|
297 | if data_conv_q is not None: |
---|
298 | output.xaxis("\\rm{Q_{x}}", output.Q_unit) |
---|
299 | output.yaxis("\\rm{Q_{y}}", output.Q_unit) |
---|
300 | else: |
---|
301 | output.xaxis("\\rm{Q_{x}}", 'A^{-1}') |
---|
302 | output.yaxis("\\rm{Q_{y}}", 'A^{-1}') |
---|
303 | if data_conv_i is not None: |
---|
304 | output.zaxis("\\rm{Intensity}", output.I_unit) |
---|
305 | else: |
---|
306 | output.zaxis("\\rm{Intensity}","cm^{-1}") |
---|
307 | |
---|
308 | # Store loading process information |
---|
309 | output.meta_data['loader'] = self.type_name |
---|
310 | |
---|
311 | return output |
---|
312 | |
---|
313 | if __name__ == "__main__": |
---|
314 | reader = Reader() |
---|
315 | print reader.read("../test/exp18_14_igor_2dqxqy.dat") |
---|
316 | |
---|