source: trunk/imports/G2pwd_fxye.py @ 1829

Last change on this file since 1829 was 1829, checked in by vondreele, 8 years ago

fix importers for HIPD & HIPPO data
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 16.8 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2015-04-30 21:34:19 +0000 (Thu, 30 Apr 2015) $
4# $Author: vondreele $
5# $Revision: 1829 $
6# $URL: trunk/imports/G2pwd_fxye.py $
7# $Id: G2pwd_fxye.py 1829 2015-04-30 21:34:19Z vondreele $
8########### SVN repository information ###################
9'''
10*Module G2pwd_fxye: GSAS data files*
11------------------------------------
12Routine to read in powder data in a variety of formats
13that are defined for GSAS.
14
15'''
16import sys
17import os.path as ospath
18import numpy as np
19import GSASIIIO as G2IO
20import GSASIIpath
21GSASIIpath.SetVersionNumber("$Revision: 1829 $")
22
23class GSAS_ReaderClass(G2IO.ImportPowderData):
24    'Routines to import powder data from a GSAS files'
25    def __init__(self):
26        super(self.__class__,self).__init__( # fancy way to self-reference
27            extensionlist=('.fxye','.raw','.gsas','.gda','gsa.','.RAW','.GSAS','.GDA','.GSA'),
28            strictExtension=False,
29            formatName = 'GSAS powder data',
30            longFormatName = 'GSAS powder data files (.fxye, .raw, .gsas...)'
31            )
32        self.clockWd = {}
33        self.TimeMap = {}
34
35    # Validate the contents -- look for a bank line
36    def ContentsValidator(self, filepointer):
37        'Validate by checking to see if the file has BANK lines'
38        #print 'ContentsValidator: '+self.formatName
39        for i,line in enumerate(filepointer):
40            self.GSAS = True
41            if i==0: # first line is always a comment
42                continue
43            if i==1 and line[:4].lower() == 'inst':
44                # 2nd line is optional instrument parameter file
45                continue
46            if line[0] == '#': continue
47            if line[:4] == 'BANK':
48                return True
49            elif line[:7] == 'Monitor': continue
50            elif line [:8] == 'TIME_MAP':          #LANSCE TOF data
51                return True
52            else:
53                self.errors = 'Unexpected information in line: '+str(i+1)
54                self.errors += '  '+str(line)
55                return False
56        return False # no bank records
57
58    def Reader(self,filename,filepointer, ParentFrame=None, **kwarg):
59        '''Read a GSAS (old formats) file of type FXY, FXYE, ESD or STD types.
60        If multiple datasets are requested, use self.repeat and buffer caching.
61        '''
62        def GetFXYEdata(File,Pos,Bank):
63            File.seek(Pos)
64            x = []
65            y = []
66            w = []
67            S = File.readline()
68            while S and S[:4] != 'BANK' and S[0] != '#':
69                vals = S.split()
70                x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
71                f = float(vals[1])
72                s = float(vals[2])
73                if f <= 0.0 or s <= 0.0:
74                    y.append(0.0)
75                    w.append(0.0)
76                else:
77                    y.append(float(vals[1]))
78                    w.append(1.0/float(vals[2])**2)
79                S = File.readline()
80            N = len(x)
81            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]   
82           
83        def GetFXYdata(File,Pos,Bank):
84            File.seek(Pos)
85            x = []
86            y = []
87            w = []
88            S = File.readline()
89            while S and S[:4] != 'BANK' and S[0] != '#':
90                vals = S.split()
91                x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
92                f = float(vals[1])
93                if f > 0.0:
94                    y.append(f)
95                    w.append(1.0/f)
96                else:             
97                    y.append(0.0)
98                    w.append(0.0)
99                S = File.readline()
100            N = len(x)
101            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
102           
103        def GetESDdata(File,Pos,Bank):
104            File.seek(Pos)
105            cons = Bank.split()
106            if 'TIME_MAP' == cons[4]:
107                start = 0
108                step = 1
109            else:
110                start = float(cons[5])/100.0               #CW: from centidegrees to degrees
111                step = float(cons[6])/100.0
112            x = []
113            y = []
114            w = []
115            S = File.readline()
116            j = 0
117            while S and S[:4] != 'BANK' and S[0] != '#':
118                if 'TIME_MAP' in S:
119                    break
120                for i in range(0,80,16):
121                    if S[i:i+8] == 8*' ':
122                        break
123                    xi = start+step*j
124                    yi = sfloat(S[i:i+8])
125                    ei = sfloat(S[i+8:i+16])
126                    x.append(xi)
127                    if yi > 0.0:
128                        y.append(yi)
129                        w.append(1.0/ei**2)
130                    else:             
131                        y.append(0.0)
132                        w.append(0.0)
133                    j += 1
134                S = File.readline()
135            N = len(x)
136            if self.clockWd:
137                x = Tmap2TOF(self.TimeMap[cons[5]],self.clockWd[cons[5]])
138            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
139       
140        def GetSTDdata(File,Pos,Bank):
141            File.seek(Pos)
142            cons = Bank.split()
143            Nch = int(cons[2])
144            if self.clockWd:
145                start = 0
146                step = 1
147            else:
148                start = float(cons[5])/100.0               #CW: from centidegrees to degrees
149                step = float(cons[6])/100.0                 #NB TOF 0.1*ms!
150            x = []
151            y = []
152            w = []
153            S = File.readline()
154            j = 0
155            while S and S[:4] != 'BANK' and S[0] != '#':
156                for i in range(0,80,8):
157                    if S[i:i+10] == 10*' ':
158                        break
159                    xi = start+step*j
160                    ni = max(sint(S[i:i+2]),1)
161                    yi = max(sfloat(S[i+2:i+8]),0.0)
162                    if yi:
163                        vi = yi/ni
164                    else:
165                        yi = 0.0
166                        vi = 0.0
167                    j += 1
168                    if j < Nch:
169                        x.append(xi)
170                        if vi <= 0.:
171                            y.append(0.)
172                            w.append(0.)
173                        else:
174                            y.append(yi)
175                            w.append(1.0/vi)
176                S = File.readline()
177            N = len(x)
178            if self.clockWd:
179                x = Tmap2TOF(self.TimeMap[cons[5]],self.clockWd[cons[5]])
180            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
181           
182        def GetALTdata(File,Pos,Bank):
183            File.seek(Pos)
184            cons = Bank.split()
185            x = []
186            y = []
187            w = []
188            S = File.readline()
189            j = 0
190            while S and S[:4] != 'BANK' and S[0] != '#':
191                for i in range(0,80,20):
192                    if S[i:i+8] == 9*' ':
193                        break
194                    xi = sfloat(S[i:i+9])/3200.
195                    yi = sfloat(S[i+9:i+16])/1000.
196                    ei = sfloat(S[i+16:i+21])/1000.
197                    x.append(xi)
198                    if yi > 0.0:
199                        y.append(yi)
200                        w.append(1.0/ei**2)
201                    else:             
202                        y.append(0.0)
203                        w.append(0.0)
204                    j += 1
205                S = File.readline()
206            N = len(x)
207            if self.clockWd:
208                x = Tmap2TOF(self.TimeMap[cons[5]],self.clockWd[cons[5]])
209            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
210           
211        def GetTimeMap(File,Pos,TimeMap):
212            File.seek(Pos)
213            cons = TimeMap[8:].split()
214            mapNo = cons[0]
215            if mapNo == '10': mapNo = '1'   #HIPD cluge!
216            Nch = int(cons[1])
217            Nrec = int(cons[2])
218            clockWd = float(cons[4])/1000.          #in mus
219            TMap = np.zeros(Nch+2,dtype=int)
220            ind = 0
221            for i in range(Nrec):
222                S = File.readline().rstrip('\n')
223                vals = S.split()
224                for val in vals:
225                    TMap[ind] = int(val)
226                    ind += 1
227            TMap = np.reshape(TMap,(-1,3))
228            TMax = TMap[-1][0]
229            Nch = TMap[-2][0]+(TMax-TMap[-2][1]+TMap[-2][2]-1)/TMap[-2][2]
230            TMap[-1] = [Nch+1,TMax,0]
231            TMap = TMap.T
232            TMap[0] -= 1
233            return TMap.T,clockWd,mapNo
234           
235        def Tmap2TOF(TMap,clockWd):
236            TOF = []
237            chWdt = []
238            Tch,T,Step = TMap[0]
239            for tmap in TMap[1:]:
240                tch,t,step = tmap
241                TOF += [T+Step*(i-Tch) for i in range(Tch,tch)]
242                Tch,T,Step = tmap
243            TOF = np.array(TOF)*clockWd
244            return TOF
245
246        x = []
247        y = []
248        w = []
249        Banks = []
250        Pos = []
251        rdbuffer = kwarg.get('buffer')
252        title = ''
253        comments = None
254#        selections = None
255
256        # reload previously saved values
257        if self.repeat and rdbuffer is not None:
258            Banks = rdbuffer.get('Banks')
259            Pos = rdbuffer.get('Pos')
260            self.selections = rdbuffer.get('selections')
261            comments = rdbuffer.get('comments')
262
263        # read through the file and find the beginning of each bank
264        # Save the offset (Pos), BANK line (Banks), comments for each bank
265        #
266        # This is going to need a fair amount of work to track line numbers
267        # in the input file.
268        if len(Banks) != len(Pos) or len(Banks) == 0:
269            try:
270                i = -1
271                while True:
272                    i += 1
273                    S = filepointer.readline()
274                    if len(S) == 0: break
275                       
276                    if i==0: # first line is always a comment
277                        self.errors = 'Error reading title'
278                        title = S[:-1]
279                        comments = [[title,]]
280                        continue
281                    if i==1 and S[:4].lower() == 'inst' and ':' in S:
282                        # 2nd line is instrument parameter file (optional)
283                        self.errors = 'Error reading instrument parameter filename'
284                        self.instparm = S.split(':')[1].strip('[]').strip()
285                        continue
286                    if S[0] == '#': # allow comments anywhere in the file
287                        # comments in fact should only preceed BANK lines
288                        comments[-1].append(S[:-1])
289                        continue
290                    if S[:4] == 'BANK':
291                        self.errors = 'Error reading bank:'
292                        self.errors += '  '+str(S)
293                        comments.append([title,])
294                        Banks.append(S)
295                        Pos.append(filepointer.tell())
296                    if S[:8] == 'TIME_MAP':     #assumes one time map; HIPPO has multiple time maps
297                        if len(Banks) == 0:
298                            self.errors = 'Error reading time map before any bank lines'
299                        else:
300                            self.errors = 'Error reading time map after bank:\n  '+str(Banks[-1])
301                        timemap,clockwd,mapNo = GetTimeMap(filepointer,filepointer.tell(),S)
302                        self.TimeMap[mapNo] = timemap
303                        self.clockWd[mapNo] = clockwd
304                       
305            except Exception as detail:
306                self.errors += '\n  '+str(detail)
307                print self.formatName+' scan error:'+str(detail) # for testing
308                import traceback
309                traceback.print_exc(file=sys.stdout)
310                return False
311
312        # Now select the bank to read
313        if not Banks: # use of ContentsValidator should prevent this error
314            print self.formatName+' scan error: no BANK records'
315            selblk = None # no block to choose
316            self.errors = 'No BANK records found (strange!)'
317            return False
318        elif len(Banks) == 1: # only one Bank, don't ask
319            selblk = 0
320        elif self.repeat and self.selections is not None:
321            # we were called to repeat the read
322            #print 'debug: repeat #',self.repeatcount,'selection',self.selections[self.repeatcount]
323            selblk = self.selections[self.repeatcount]
324            self.repeatcount += 1
325            if self.repeatcount >= len(self.selections): self.repeat = False
326        else:                       # choose from options
327            if not len(self.selections):    #use previous selection, otherwise...
328                self.selections = self.MultipleBlockSelector(
329                    Banks,
330                    ParentFrame=ParentFrame,
331                    title='Select Bank(s) to read from the list below',
332                    size=(600,100),
333                    header='Dataset Selector')
334            if len(self.selections) == 0: return False
335            selblk = self.selections[0] # select first in list
336            if len(self.selections) > 1: # prepare to loop through again
337                self.repeat = True
338                self.repeatcount = 1
339                if rdbuffer is not None:
340                    rdbuffer['Banks'] = Banks
341                    rdbuffer['Pos'] = Pos
342                    rdbuffer['selections'] = self.selections
343                    rdbuffer['comments'] = comments
344
345        # got a selection, now read it
346        Bank = Banks[selblk]
347        try:
348            bnkNo = int(Bank.split()[1])
349        except ValueError:
350            bnkNo = 1
351        try:
352            if 'FXYE' in Bank:
353                self.errors = 'Error reading FXYE data in Bank\n  '+Banks[selblk]
354                self.powderdata = GetFXYEdata(filepointer,Pos[selblk],Banks[selblk])
355            elif 'FXY' in Bank:
356                self.errors = 'Error reading FXY data in Bank\n  '+Banks[selblk]
357                self.powderdata = GetFXYdata(filepointer,Pos[selblk],Banks[selblk])
358            elif 'ESD' in Bank:
359                self.errors = 'Error reading ESD data in Bank\n  '+Banks[selblk]
360                self.powderdata = GetESDdata(filepointer,Pos[selblk],Banks[selblk])
361            elif 'STD' in Bank:
362                self.errors = 'Error reading STD data in Bank\n  '+Banks[selblk]
363                self.powderdata = GetSTDdata(filepointer,Pos[selblk],Banks[selblk])
364            elif 'ALT' in Bank:
365                self.errors = 'Error reading ALT data in Bank\n  '+Banks[selblk]
366                self.powderdata = GetALTdata(filepointer,Pos[selblk],Banks[selblk])
367            else:
368                self.errors = 'Error reading STD data in Bank\n  '+Banks[selblk]
369                self.powderdata = GetSTDdata(filepointer,Pos[selblk],Banks[selblk])
370        except Exception as detail:
371            self.errors += '\n  '+str(detail)
372            print self.formatName+' read error:'+str(detail) # for testing
373            import traceback
374            traceback.print_exc(file=sys.stdout)
375            return False
376
377        self.errors = 'Error processing information after read complete'
378        if comments is not None:
379            self.comments = comments[selblk]
380        self.powderentry[0] = filename
381        self.powderentry[1] = Pos # position offset (never used, I hope)
382        self.powderentry[2] = bnkNo #selblk+1 # bank number
383        self.idstring = ospath.basename(filename) + ' Bank '+str(bnkNo) #selblk+1)
384        self.numbanks=len(Banks)
385        # scan comments for temperature & radius
386        Temperature = 300.
387        for S in self.comments:
388            if 'Temp' in S.split('=')[0]:
389                try:
390                    Temperature = float(S.split('=')[1])
391                except:
392                    pass
393            elif 'Gonio' in S.split('=')[0]:
394                try:
395                    self.Sample['Gonio. radius'] = float(S.split('=')[1])
396                except:
397                    pass
398            elif 'Omega' in S.split('=')[0] or 'Theta' in S.split('=')[0]:  #HIPD weirdness
399                try:
400                    self.Sample['Omega'] = 90.-float(S.split('=')[1])
401                except:
402                    pass
403            elif 'Chi' in S.split('=')[0]:
404                try:
405                    self.Sample['Chi'] = -float(S.split('=')[1])
406                except:
407                    pass                   
408            elif 'Phi' in S.split('=')[0]:
409                try:
410                    self.Sample['Phi'] = float(S.split('=')[1])
411                except:
412                    pass                   
413        self.Sample['Temperature'] = Temperature
414        return True       
415
416def sfloat(S):
417    'convert a string to a float, treating an all-blank string as zero'
418    if S.strip():
419        return float(S)
420    else:
421        return 0.0
422
423def sint(S):
424    'convert a string to an integer, treating an all-blank string as zero'
425    if S.strip():
426        return int(S)
427    else:
428        return 0
Note: See TracBrowser for help on using the repository browser.