source: trunk/imports/G2pwd_CIF.py @ 1074

Last change on this file since 1074 was 1074, checked in by toby, 8 years ago

more CIF work

  • Property svn:eol-style set to native
File size: 13.6 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2012-02-13 11:33:35 -0600 (Mon, 13 Feb 2012) $
4# $Author: vondreele & toby $
5# $Revision: 482 $
6# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/G2importphase.py $
7# $Id: G2importphase.py 482 2012-02-13 17:33:35Z vondreele $
8########### SVN repository information ###################
9# routines to read in structure factors from a CIF
10#
11import sys
12import numpy as np
13import os.path
14import GSASIIIO as G2IO
15import CifFile as cif # PyCifRW from James Hester
16import GSASIIpath
17GSASIIpath.SetVersionNumber("$Revision: 810 $")
18
19class CIFpwdReader(G2IO.ImportPowderData):
20    'Routines to import powder data from a CIF file'
21    def __init__(self):
22        super(self.__class__,self).__init__( # fancy way to self-reference
23            extensionlist=('.CIF','.cif'),
24            strictExtension=False,
25            formatName = 'CIF',
26            longFormatName = 'Powder data from CIF'
27            )
28    # Validate the contents
29    def ContentsValidator(self, filepointer):
30        for i,line in enumerate(filepointer):
31            if i >= 1000: break
32            ''' Encountered only blank lines or comments in first 1000
33            lines. This is unlikely, but assume it is CIF since we are
34            even less likely to find a file with nothing but hashes and
35            blank lines'''
36            line = line.strip()
37            if len(line) == 0:
38                continue # ignore blank lines
39            elif line.startswith('#'):
40                continue # ignore comments
41            elif line.startswith('data_'):
42                return True
43            else:
44                return False # found something else
45        return True
46
47    def Reader(self,filename,filepointer, ParentFrame=None, **kwarg):
48        # Define lists of data names used for holding powder diffraction data
49        xDataItems = (  # "x-axis" data names"
50            ('_pd_meas_2theta_range_min', '_pd_meas_2theta_range_max', '_pd_meas_2theta_range_inc'),
51            ('_pd_proc_2theta_range_min', '_pd_proc_2theta_range_max', '_pd_proc_2theta_range_inc'),
52            '_pd_meas_2theta_scan',
53            '_pd_meas_time_of_flight',
54            '_pd_proc_2theta_corrected',
55            #'_pd_proc_d_spacing',
56            #'_pd_proc_energy_incident',
57            #'_pd_proc_energy_detection',
58            #'_pd_proc_recip_len_q',
59            #'_pd_proc_wavelength',
60        )
61        intDataItems = ( # intensity axis data names
62            '_pd_meas_counts_total',
63            #'_pd_meas_counts_background',
64            #'_pd_meas_counts_container',
65            '_pd_meas_intensity_total',
66            #'_pd_meas_intensity_background',
67            #'_pd_meas_intensity_container',
68            '_pd_proc_intensity_net',
69            '_pd_proc_intensity_total',
70            #'_pd_proc_intensity_bkg_calc',
71            #'_pd_proc_intensity_bkg_fix',
72            '_pd_calc_intensity_net',   # allow computed patterns as input data?
73            '_pd_calc_intensity_total',
74            )
75
76        ESDDataItems = ( # items that can be used to compute uncertainties for obs data
77            '_pd_proc_ls_weight',
78            '_pd_meas_counts_total'
79            )
80       
81        ModDataItems = ( # items that modify the use of the data
82            '_pd_meas_step_count_time',
83            '_pd_meas_counts_monitor',
84            '_pd_meas_intensity_monitor',
85            '_pd_proc_intensity_norm',
86            '_pd_proc_intensity_incident',
87            )
88        rdbuffer = kwarg.get('buffer')
89        cf = None
90        choicelist = None
91        selections = None
92        # reload previously saved values
93        if self.repeat and rdbuffer is not None:
94            cf = rdbuffer.get('lastcif')
95            choicelist = rdbuffer.get('choicelist')
96            print 'debug: Reuse previously parsed CIF'
97            selections = rdbuffer.get('selections')
98        try:
99            if cf is None:
100                self.ShowBusy() # this can take a while
101                cf = G2IO.ReadCIF(filename)
102                self.DoneBusy()
103        except Exception as detail:
104            print self.formatName+' read error:'+str(detail) # for testing
105            import traceback
106            traceback.print_exc(file=sys.stdout)
107            return False
108               
109        # scan all blocks for sets of data
110        if choicelist is None:
111            choicelist = []
112            for blk in cf.keys():
113                blkkeys = [k.lower() for k in cf[blk].keys()] # save a list of the data items, since we will use it often
114                # scan through block for x items
115                xldict = {}
116                for x in xDataItems:
117                    if type(x) is tuple: # check for the presence of all three items that define a range of data
118                        if not all([i in blkkeys for i in x]): continue
119                        try:
120                            items = [float(cf[blk][xi]) for xi in x]
121                            l = 1 + int(0.5 + (items[1]-items[0])/items[2])
122                        except:
123                            continue
124                    else:
125                        if x not in blkkeys: continue
126                        l = len(cf[blk][x])
127                    if xldict.get(l) is None:
128                        xldict[l] = [x]
129                    else:
130                        xldict[l].append(x)
131                # now look for matching intensity items
132                yldict = {}
133                suldict = {}
134                for y in intDataItems:
135                    if y in blkkeys:
136                        l = len(cf[blk][y])
137                        if yldict.get(l) is None:
138                            yldict[l] = [y]
139                        else:
140                            yldict[l].append(y)
141                        # now check if the first item has an uncertainty
142                        if cif.get_number_with_esd(cf[blk][y][0])[1] is None: continue
143                        if suldict.get(l) is None:
144                            suldict[l] = [y]
145                        else:
146                            suldict[l].append(y)
147                for y in ESDDataItems:
148                    if y in blkkeys:
149                        l = len(cf[blk][y])
150                        if suldict.get(l) is None:
151                            suldict[l] = [y]
152                        else:
153                            suldict[l].append(y)
154                modldict = {}
155                for y in ModDataItems:
156                    if y in blkkeys:
157                        l = len(cf[blk][y])
158                        if modldict.get(l) is None:
159                            modldict[l] = [y]
160                        else:
161                            modldict[l].append(y)
162                for l in xldict:
163                    if yldict.get(l) is None: continue
164                    choicelist.append([blk,l,xldict[l],yldict[l],suldict.get(l,[]),modldict.get(l,[])])
165                    #print blk,l,xldict[l],yldict[l],suldict.get(l,[]),modldict.get(l,[])
166        if not choicelist:
167            selblk = None # no block to choose
168            return False
169        elif len(choicelist) == 1: # only one choice
170            selblk = 0
171        elif self.repeat and selections is not None:
172            # we were called to repeat the read
173            print 'debug: repeat #',self.repeatcount,'selection',selections[self.repeatcount]
174            selblk = selections[self.repeatcount]
175            self.repeatcount += 1
176            if self.repeatcount >= len(selections): self.repeat = False
177        else:                       # choose from options
178            # compile a list of choices for the user
179            choices = []
180            for blk,l,x,y,su,mod in choicelist:
181                sx = x[0]
182                if len(x) > 1: sx += '...'
183                sy = y[0]
184                if len(y) > 1: sy += '...'
185                choices.append(
186                    'Block '+str(blk)+', '+str(l)+' points. X='+sx+' & Y='+sy
187                    )
188            selections = self.MultipleBlockSelector(
189                choices,
190                ParentFrame=ParentFrame,
191                title='Select dataset(s) to read from the list below',
192                size=(600,100),
193                header='Dataset Selector')
194            if len(selections) == 0: return False
195            selblk = selections[0] # select first in list
196            if len(selections) > 1: # prepare to loop through again
197                self.repeat = True
198                self.repeatcount = 1
199                if rdbuffer is not None:
200                    rdbuffer['selections'] = selections
201                    rdbuffer['lastcif'] = cf # save the parsed cif for the next loop
202                    rdbuffer['choicelist'] = choicelist # save the parsed choices for the future
203
204        # got a selection, now read it
205        blk,l,xch,ych,such,modch = choicelist[selblk]
206        xi,yi,sui,modi = 0,0,0,0
207        if len(xch) > 1 or len(ych) > 1 or len(such) > 1 or len(modch) > 0:
208            choices = []
209            chlbls = []
210            chlbls.append('Select the scanned data item')
211            xchlst = []
212            for x in xch:
213                if type(x) is tuple:
214                    xchlst.append(x[0])
215                else:
216                    xchlst.append(x)
217            choices.append(xchlst)
218            chlbls.append('Select the intensity data item')
219            choices.append(ych)
220            chlbls.append('Select the data item to be used for weighting')
221            choices.append(such)
222            chlbls.append('Divide intensities by data item')
223            choices.append(['none']+modch)
224            res = self.MultipleChoicesDialog(choices,chlbls)
225            if not res: return False
226            xi,yi,sui,modi = res
227
228        # now read in the values
229        # x-values
230        xcf = xch[xi]
231        if type(xcf) is tuple:
232            vals = [float(cf[blk][xi]) for xi in xcf]
233            x = np.array(
234                [(i*vals[2] + vals[0]) for i in range(1 + int(0.5 + (vals[1]-vals[0])/vals[2]))]
235                )
236        else:
237            vl = []
238            for val in cf[blk].get(xcf,'?'):
239                v,e = cif.get_number_with_esd(val)
240                if v is None: # not parsed
241                    vl.append(np.NaN)
242                else:
243                    vl.append(v)
244            x = np.array(vl)
245        # y-values
246        ycf = ych[yi]
247        vl = []
248        for val in cf[blk].get(ycf,'?'):
249            v,e = cif.get_number_with_esd(val)
250            if v is None: # not parsed
251                vl.append(np.NaN)
252            else:
253                vl.append(v)
254        y = np.array(vl)
255        # weights
256        if sui == -1:
257            # no weights
258            vl = np.zeros(len(x)) + 1.
259        else:
260            sucf = such[sui]
261            if sucf ==  '_pd_proc_ls_weight':
262                for val in cf[blk].get(sucf,'?'):
263                    v,e = cif.get_number_with_esd(val)
264                    if v is None: # not parsed
265                        vl.append(0.)
266                    else:
267                        vl.append(v)
268            elif sucf ==  '_pd_meas_counts_total':
269                for val in cf[blk].get(sucf,'?'):
270                    v,e = cif.get_number_with_esd(val)
271                    if v is None: # not parsed
272                        vl.append(0.)
273                    elif v <= 0:
274                        vl.append(1.)
275                    else:
276                        vl.append(1./v)
277            else:
278                for val in cf[blk].get(sucf,'?'):
279                    v,e = cif.get_number_with_esd(val)
280                    if v is None: # not parsed
281                        vl.append(0.0)
282                    elif v <= 0:
283                        vl.append(0.0)
284                    else:
285                        vl.append(1./(v*v))
286        w = np.array(vl)
287        # intensity modification factor
288        if modi >= 1:
289            ycf = modch[modi-1]
290            vl = []
291            for val in cf[blk].get(ycf,'?'):
292                v,e = cif.get_number_with_esd(val)
293                if v is None: # not parsed
294                    vl.append(np.NaN)
295                else:
296                    vl.append(v)
297            y /= np.array(vl)
298            w /= np.array(vl)
299        N = len(x)
300        self.powderdata = [
301                np.array(x), # x-axis values
302                np.array(y), # powder pattern intensities
303                np.array(w), # 1/sig(intensity)^2 values (weights)
304                np.zeros(N), # calc. intensities (zero)
305                np.zeros(N), # calc. background (zero)
306                np.zeros(N), # obs-calc profiles
307            ]
308        self.powderentry[0] = filename
309        #self.powderentry[1] = pos # bank offset (N/A here)
310        self.powderentry[2] = 1 # xye file only has one bank
311        self.idstring = os.path.basename(filename) + ': ' + blk
312        if cf[blk].get('_diffrn_radiation_probe'):
313            if cf[blk]['_diffrn_radiation_probe'] == 'neutron':
314                self.instdict['type'] = 'SNC'
315                #if cf[blk].get('_pd_meas_time_of_flight'): self.instdict['type'] = 'SNT' # not supported yet
316            else:
317                self.instdict['type'] = 'SXC'
318        if cf[blk].get('_diffrn_radiation_wavelength'):
319            val = cf[blk]['_diffrn_radiation_wavelength']
320            wl = []
321            if type(val) is list:
322                for v in val:
323                    w,e = cif.get_number_with_esd(v)
324                    if w: wl.append(w)
325            else:
326                w,e = cif.get_number_with_esd(val)
327                if w: wl.append(w)
328            if wl: self.instdict['wave'] = wl
329        if cf[blk].get('_diffrn_ambient_temperature'):
330            val = cf[blk]['_diffrn_ambient_temperature']
331            w,e = cif.get_number_with_esd(val)
332            if w: 
333                self.Sample['Temperature'] = w
334        return True
Note: See TracBrowser for help on using the repository browser.