[959eb01] | 1 | """ |
---|
| 2 | DANSE/SANS file reader |
---|
| 3 | """ |
---|
| 4 | ############################################################################ |
---|
| 5 | #This software was developed by the University of Tennessee as part of the |
---|
| 6 | #Distributed Data Analysis of Neutron Scattering Experiments (DANSE) |
---|
| 7 | #project funded by the US National Science Foundation. |
---|
| 8 | #If you use DANSE applications to do scientific research that leads to |
---|
| 9 | #publication, we ask that you acknowledge the use of the software with the |
---|
| 10 | #following sentence: |
---|
| 11 | #This work benefited from DANSE software developed under NSF award DMR-0520547. |
---|
| 12 | #copyright 2008, University of Tennessee |
---|
| 13 | ############################################################################# |
---|
| 14 | import math |
---|
| 15 | import os |
---|
| 16 | import sys |
---|
| 17 | import numpy as np |
---|
| 18 | import logging |
---|
| 19 | from sas.sascalc.dataloader.data_info import Data2D, Detector |
---|
| 20 | from sas.sascalc.dataloader.manipulations import reader2D_converter |
---|
| 21 | |
---|
| 22 | logger = logging.getLogger(__name__) |
---|
| 23 | |
---|
| 24 | # Look for unit converter |
---|
| 25 | has_converter = True |
---|
| 26 | try: |
---|
| 27 | from sas.sascalc.data_util.nxsunit import Converter |
---|
| 28 | except: |
---|
| 29 | has_converter = False |
---|
| 30 | |
---|
| 31 | |
---|
| 32 | class Reader: |
---|
| 33 | """ |
---|
| 34 | Example data manipulation |
---|
| 35 | """ |
---|
| 36 | ## File type |
---|
| 37 | type_name = "DANSE" |
---|
| 38 | ## Wildcards |
---|
| 39 | type = ["DANSE files (*.sans)|*.sans"] |
---|
| 40 | ## Extension |
---|
| 41 | ext = ['.sans', '.SANS'] |
---|
| 42 | |
---|
| 43 | def read(self, filename=None): |
---|
| 44 | """ |
---|
| 45 | Open and read the data in a file |
---|
| 46 | @param file: path of the file |
---|
| 47 | """ |
---|
| 48 | |
---|
| 49 | read_it = False |
---|
| 50 | for item in self.ext: |
---|
| 51 | if filename.lower().find(item) >= 0: |
---|
| 52 | read_it = True |
---|
| 53 | |
---|
| 54 | if read_it: |
---|
| 55 | try: |
---|
| 56 | datafile = open(filename, 'r') |
---|
| 57 | except: |
---|
| 58 | raise RuntimeError,"danse_reader cannot open %s" % (filename) |
---|
| 59 | |
---|
| 60 | # defaults |
---|
| 61 | # wavelength in Angstrom |
---|
| 62 | wavelength = 10.0 |
---|
| 63 | # Distance in meter |
---|
| 64 | distance = 11.0 |
---|
| 65 | # Pixel number of center in x |
---|
| 66 | center_x = 65 |
---|
| 67 | # Pixel number of center in y |
---|
| 68 | center_y = 65 |
---|
| 69 | # Pixel size [mm] |
---|
| 70 | pixel = 5.0 |
---|
| 71 | # Size in x, in pixels |
---|
| 72 | size_x = 128 |
---|
| 73 | # Size in y, in pixels |
---|
| 74 | size_y = 128 |
---|
| 75 | # Format version |
---|
| 76 | fversion = 1.0 |
---|
| 77 | |
---|
| 78 | output = Data2D() |
---|
| 79 | output.filename = os.path.basename(filename) |
---|
| 80 | detector = Detector() |
---|
| 81 | output.detector.append(detector) |
---|
| 82 | |
---|
| 83 | output.data = np.zeros([size_x,size_y]) |
---|
| 84 | output.err_data = np.zeros([size_x, size_y]) |
---|
| 85 | |
---|
| 86 | data_conv_q = None |
---|
| 87 | data_conv_i = None |
---|
| 88 | |
---|
| 89 | if has_converter == True and output.Q_unit != '1/A': |
---|
| 90 | data_conv_q = Converter('1/A') |
---|
| 91 | # Test it |
---|
| 92 | data_conv_q(1.0, output.Q_unit) |
---|
| 93 | |
---|
| 94 | if has_converter == True and output.I_unit != '1/cm': |
---|
| 95 | data_conv_i = Converter('1/cm') |
---|
| 96 | # Test it |
---|
| 97 | data_conv_i(1.0, output.I_unit) |
---|
| 98 | |
---|
| 99 | read_on = True |
---|
| 100 | while read_on: |
---|
| 101 | line = datafile.readline() |
---|
| 102 | if line.find("DATA:") >= 0: |
---|
| 103 | read_on = False |
---|
| 104 | break |
---|
| 105 | toks = line.split(':') |
---|
| 106 | if toks[0] == "FORMATVERSION": |
---|
| 107 | fversion = float(toks[1]) |
---|
| 108 | if toks[0] == "WAVELENGTH": |
---|
| 109 | wavelength = float(toks[1]) |
---|
| 110 | elif toks[0] == "DISTANCE": |
---|
| 111 | distance = float(toks[1]) |
---|
| 112 | elif toks[0] == "CENTER_X": |
---|
| 113 | center_x = float(toks[1]) |
---|
| 114 | elif toks[0] == "CENTER_Y": |
---|
| 115 | center_y = float(toks[1]) |
---|
| 116 | elif toks[0] == "PIXELSIZE": |
---|
| 117 | pixel = float(toks[1]) |
---|
| 118 | elif toks[0] == "SIZE_X": |
---|
| 119 | size_x = int(toks[1]) |
---|
| 120 | elif toks[0] == "SIZE_Y": |
---|
| 121 | size_y = int(toks[1]) |
---|
| 122 | |
---|
| 123 | # Read the data |
---|
| 124 | data = [] |
---|
| 125 | error = [] |
---|
| 126 | if fversion == 1.0: |
---|
| 127 | data_str = datafile.readline() |
---|
| 128 | data = data_str.split(' ') |
---|
| 129 | else: |
---|
| 130 | read_on = True |
---|
| 131 | while read_on: |
---|
| 132 | data_str = datafile.readline() |
---|
| 133 | if len(data_str) == 0: |
---|
| 134 | read_on = False |
---|
| 135 | else: |
---|
| 136 | toks = data_str.split() |
---|
| 137 | try: |
---|
| 138 | val = float(toks[0]) |
---|
| 139 | err = float(toks[1]) |
---|
| 140 | if data_conv_i is not None: |
---|
| 141 | val = data_conv_i(val, units=output._yunit) |
---|
| 142 | err = data_conv_i(err, units=output._yunit) |
---|
| 143 | data.append(val) |
---|
| 144 | error.append(err) |
---|
| 145 | except: |
---|
| 146 | logger.info("Skipping line:%s,%s" %(data_str, |
---|
| 147 | sys.exc_value)) |
---|
| 148 | |
---|
| 149 | # Initialize |
---|
| 150 | x_vals = [] |
---|
| 151 | y_vals = [] |
---|
| 152 | ymin = None |
---|
| 153 | ymax = None |
---|
| 154 | xmin = None |
---|
| 155 | xmax = None |
---|
| 156 | |
---|
| 157 | # Qx and Qy vectors |
---|
| 158 | theta = pixel / distance / 100.0 |
---|
| 159 | stepq = 4.0 * math.pi / wavelength * math.sin(theta / 2.0) |
---|
| 160 | for i_x in range(size_x): |
---|
| 161 | theta = (i_x - center_x + 1) * pixel / distance / 100.0 |
---|
| 162 | qx = 4.0 * math.pi / wavelength * math.sin(theta / 2.0) |
---|
| 163 | |
---|
| 164 | if has_converter == True and output.Q_unit != '1/A': |
---|
| 165 | qx = data_conv_q(qx, units=output.Q_unit) |
---|
| 166 | |
---|
| 167 | x_vals.append(qx) |
---|
[235f514] | 168 | if xmin is None or qx < xmin: |
---|
[959eb01] | 169 | xmin = qx |
---|
[235f514] | 170 | if xmax is None or qx > xmax: |
---|
[959eb01] | 171 | xmax = qx |
---|
| 172 | |
---|
| 173 | ymin = None |
---|
| 174 | ymax = None |
---|
| 175 | for i_y in range(size_y): |
---|
| 176 | theta = (i_y - center_y + 1) * pixel / distance / 100.0 |
---|
| 177 | qy = 4.0 * math.pi / wavelength * math.sin(theta/2.0) |
---|
| 178 | |
---|
| 179 | if has_converter == True and output.Q_unit != '1/A': |
---|
| 180 | qy = data_conv_q(qy, units=output.Q_unit) |
---|
| 181 | |
---|
| 182 | y_vals.append(qy) |
---|
[235f514] | 183 | if ymin is None or qy < ymin: |
---|
[959eb01] | 184 | ymin = qy |
---|
[235f514] | 185 | if ymax is None or qy > ymax: |
---|
[959eb01] | 186 | ymax = qy |
---|
| 187 | |
---|
| 188 | # Store the data in the 2D array |
---|
| 189 | i_x = 0 |
---|
| 190 | i_y = -1 |
---|
| 191 | |
---|
| 192 | for i_pt in range(len(data)): |
---|
| 193 | try: |
---|
| 194 | value = float(data[i_pt]) |
---|
| 195 | except: |
---|
| 196 | # For version 1.0, the data were still |
---|
| 197 | # stored as strings at this point. |
---|
| 198 | msg = "Skipping entry (v1.0):%s,%s" % (str(data[i_pt]), |
---|
| 199 | sys.exc_value) |
---|
| 200 | logger.info(msg) |
---|
| 201 | |
---|
| 202 | # Get bin number |
---|
| 203 | if math.fmod(i_pt, size_x) == 0: |
---|
| 204 | i_x = 0 |
---|
| 205 | i_y += 1 |
---|
| 206 | else: |
---|
| 207 | i_x += 1 |
---|
| 208 | |
---|
| 209 | output.data[i_y][i_x] = value |
---|
| 210 | if fversion>1.0: |
---|
| 211 | output.err_data[i_y][i_x] = error[i_pt] |
---|
| 212 | |
---|
| 213 | # Store all data |
---|
| 214 | # Store wavelength |
---|
| 215 | if has_converter == True and output.source.wavelength_unit != 'A': |
---|
| 216 | conv = Converter('A') |
---|
| 217 | wavelength = conv(wavelength, |
---|
| 218 | units=output.source.wavelength_unit) |
---|
| 219 | output.source.wavelength = wavelength |
---|
| 220 | |
---|
| 221 | # Store distance |
---|
| 222 | if has_converter == True and detector.distance_unit != 'm': |
---|
| 223 | conv = Converter('m') |
---|
| 224 | distance = conv(distance, units=detector.distance_unit) |
---|
| 225 | detector.distance = distance |
---|
| 226 | |
---|
| 227 | # Store pixel size |
---|
| 228 | if has_converter == True and detector.pixel_size_unit != 'mm': |
---|
| 229 | conv = Converter('mm') |
---|
| 230 | pixel = conv(pixel, units=detector.pixel_size_unit) |
---|
| 231 | detector.pixel_size.x = pixel |
---|
| 232 | detector.pixel_size.y = pixel |
---|
| 233 | |
---|
| 234 | # Store beam center in distance units |
---|
| 235 | detector.beam_center.x = center_x * pixel |
---|
| 236 | detector.beam_center.y = center_y * pixel |
---|
| 237 | |
---|
| 238 | # Store limits of the image (2D array) |
---|
| 239 | xmin = xmin - stepq / 2.0 |
---|
| 240 | xmax = xmax + stepq / 2.0 |
---|
| 241 | ymin = ymin - stepq /2.0 |
---|
| 242 | ymax = ymax + stepq / 2.0 |
---|
| 243 | |
---|
| 244 | if has_converter == True and output.Q_unit != '1/A': |
---|
| 245 | xmin = data_conv_q(xmin, units=output.Q_unit) |
---|
| 246 | xmax = data_conv_q(xmax, units=output.Q_unit) |
---|
| 247 | ymin = data_conv_q(ymin, units=output.Q_unit) |
---|
| 248 | ymax = data_conv_q(ymax, units=output.Q_unit) |
---|
| 249 | output.xmin = xmin |
---|
| 250 | output.xmax = xmax |
---|
| 251 | output.ymin = ymin |
---|
| 252 | output.ymax = ymax |
---|
| 253 | |
---|
| 254 | # Store x and y axis bin centers |
---|
| 255 | output.x_bins = x_vals |
---|
| 256 | output.y_bins = y_vals |
---|
| 257 | |
---|
| 258 | # Units |
---|
| 259 | if data_conv_q is not None: |
---|
| 260 | output.xaxis("\\rm{Q_{x}}", output.Q_unit) |
---|
| 261 | output.yaxis("\\rm{Q_{y}}", output.Q_unit) |
---|
| 262 | else: |
---|
| 263 | output.xaxis("\\rm{Q_{x}}", 'A^{-1}') |
---|
| 264 | output.yaxis("\\rm{Q_{y}}", 'A^{-1}') |
---|
| 265 | |
---|
| 266 | if data_conv_i is not None: |
---|
| 267 | output.zaxis("\\rm{Intensity}", output.I_unit) |
---|
| 268 | else: |
---|
| 269 | output.zaxis("\\rm{Intensity}", "cm^{-1}") |
---|
| 270 | |
---|
| 271 | if not fversion >= 1.0: |
---|
| 272 | msg = "Danse_reader can't read this file %s" % filename |
---|
| 273 | raise ValueError, msg |
---|
| 274 | else: |
---|
| 275 | logger.info("Danse_reader Reading %s \n" % filename) |
---|
| 276 | |
---|
| 277 | # Store loading process information |
---|
| 278 | output.meta_data['loader'] = self.type_name |
---|
| 279 | output = reader2D_converter(output) |
---|
| 280 | return output |
---|
| 281 | |
---|
| 282 | return None |
---|