source: trunk/imports/G2pwd_CIF.py @ 3098

Last change on this file since 3098 was 3098, checked in by toby, 6 years ago

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