source: trunk/imports/G2pwd_fxye.py @ 1476

Last change on this file since 1476 was 1476, checked in by vondreele, 7 years ago

fix handling of zero powder profile intensities; set w=0. (not 1.)

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