Changeset 1443


Ignore:
Timestamp:
Jul 28, 2014 4:53:56 PM (7 years ago)
Author:
vondreele
Message:

add calibration of lam, difC, etc. from index peak positions
new plotCalib routine

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIgrid.py

    r1442 r1443  
    113113
    114114[ wxID_INSTPRMRESET,wxID_CHANGEWAVETYPE,wxID_INSTCOPY, wxID_INSTFLAGCOPY, wxID_INSTLOAD,
    115     wxID_INSTSAVE, wxID_INST1VAL
    116 ] = [wx.NewId() for item in range(7)]
     115    wxID_INSTSAVE, wxID_INST1VAL, wxID_INSTCALIB,
     116] = [wx.NewId() for item in range(8)]
    117117
    118118[ wxID_UNDO,wxID_LSQPEAKFIT,wxID_LSQONECYCLE,wxID_RESETSIGGAM,wxID_CLEARPEAKS,wxID_AUTOSEARCH,
     
    28712871        self.InstEdit = wx.Menu(title='')
    28722872        self.InstMenu.Append(menu=self.InstEdit, title='Operations')
     2873        self.InstEdit.Append(help='Calibrate from indexed peaks',
     2874            id=wxID_INSTCALIB, kind=wx.ITEM_NORMAL,text='Calibrate')
    28732875        self.InstEdit.Append(help='Reset instrument profile parameters to default',
    28742876            id=wxID_INSTLOAD, kind=wx.ITEM_NORMAL,text='Load profile...')
  • trunk/GSASIIlattice.py

    r1439 r1443  
    428428def Dsp2pos(Inst,dsp):
    429429    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
    430     ignores secondary effects (e.g. difA in TOF) - maybe later?
    431430    '''
    432431    if 'C' in Inst['Type'][0]:
     
    435434    else:   #'T'OF
    436435        pos = Inst['difC'][1]*dsp+Inst['Zero'][1]+Inst['difA'][1]*dsp**2+Inst.get('difB',[0,0,False])[1]*dsp**3
    437     return pos             
     436    return pos
     437   
     438def getPeakPos(dataType,parmdict,dsp):
     439    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
     440    '''
     441    if 'C' in dataType:
     442        pos = 2.0*asind(parmdict['Lam']/(2.*dsp))+parmdict['Zero']
     443    else:   #'T'OF
     444        pos = parmdict['difC']*dsp+parmdict['difA']*dsp**2+parmdict['difB']*dsp**3+parmdict['Zero']
     445    return pos
     446                   
    438447   
    439448def calc_rDsq(H,A):
  • trunk/GSASIIplot.py

    r1442 r1443  
    18281828       
    18291829################################################################################
     1830##### PlotCalib
     1831################################################################################
     1832           
     1833def PlotCalib(G2frame,Inst,XY,newPlot=False):
     1834    '''plot of CW or TOF peak calibration
     1835    '''
     1836    def OnMotion(event):
     1837        xpos = event.xdata
     1838        if xpos:                                        #avoid out of frame mouse position
     1839            ypos = event.ydata
     1840            Page.canvas.SetCursor(wx.CROSS_CURSOR)
     1841            try:
     1842                G2frame.G2plotNB.status.SetStatusText('X =%9.3f %s =%9.3f'%(xpos,Title,ypos),1)                   
     1843            except TypeError:
     1844                G2frame.G2plotNB.status.SetStatusText('Select '+Title+' pattern first',1)
     1845
     1846    Title = 'Position calibration'
     1847    try:
     1848        plotNum = G2frame.G2plotNB.plotList.index(Title)
     1849        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
     1850        if not newPlot:
     1851            Plot = Page.figure.gca()
     1852            xylim = Plot.get_xlim(),Plot.get_ylim()
     1853        Page.figure.clf()
     1854        Plot = Page.figure.gca()
     1855    except ValueError:
     1856        newPlot = True
     1857        Plot = G2frame.G2plotNB.addMpl(Title).gca()
     1858        plotNum = G2frame.G2plotNB.plotList.index(Title)
     1859        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
     1860        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
     1861   
     1862    Page.Choice = None
     1863    Page.SetFocus()
     1864    G2frame.G2plotNB.status.DestroyChildren()
     1865    Plot.set_title(Title)
     1866    Plot.set_xlabel(r'd-spacing',fontsize=14)
     1867    if 'C' in Inst['Type'][0]:
     1868        Plot.set_ylabel(r'$\mathsf{\Delta(2\theta)}$',fontsize=14)
     1869    else:
     1870        Plot.set_ylabel(r'$\mathsf{\Delta}T/T$',fontsize=14)
     1871    colors=['b','g','r','c','m','k']
     1872    for ixy,xy in enumerate(XY):
     1873        X,Y = xy
     1874        Yc = G2lat.Dsp2pos(Inst,X)
     1875        if 'C' in Inst['Type'][0]:
     1876            Y = Y-Yc
     1877        else:
     1878            Y = (Y-Yc)/Yc
     1879        Plot.plot(X,Y,'k+',picker=False)
     1880    if not newPlot:
     1881        Page.toolbar.push_current()
     1882        Plot.set_xlim(xylim[0])
     1883        Plot.set_ylim(xylim[1])
     1884        xylim = []
     1885        Page.toolbar.push_current()
     1886        Page.toolbar.draw()
     1887    else:
     1888        Page.canvas.draw()
     1889
     1890################################################################################
    18301891##### PlotXY
    18311892################################################################################
  • trunk/GSASIIpwd.py

    r1439 r1443  
    11231123    backVary += peaksVary   
    11241124    return bakType,backDict,backVary
     1125   
     1126def DoCalibInst(IndexPeaks,Inst):
     1127   
     1128    def SetInstParms():
     1129        dataType = Inst['Type'][0]
     1130        insVary = []
     1131        insNames = []
     1132        insVals = []
     1133        for parm in Inst:
     1134            insNames.append(parm)
     1135            insVals.append(Inst[parm][1])
     1136            if parm in ['Lam','difC','difA','difB','Zero',]:
     1137                if Inst[parm][2]:
     1138                    insVary.append(parm)
     1139        instDict = dict(zip(insNames,insVals))
     1140        return dataType,instDict,insVary
     1141       
     1142    def GetInstParms(parmDict,Inst,varyList):
     1143        for name in Inst:
     1144            Inst[name][1] = parmDict[name]
     1145       
     1146    def InstPrint(Inst,sigDict):
     1147        print 'Instrument Parameters:'
     1148        if 'C' in Inst['Type'][0]:
     1149            ptfmt = "%12.6f"
     1150        else:
     1151            ptfmt = "%12.3f"
     1152        ptlbls = 'names :'
     1153        ptstr =  'values:'
     1154        sigstr = 'esds  :'
     1155        for parm in Inst:
     1156            if parm in  ['Lam','difC','difA','difB','Zero',]:
     1157                ptlbls += "%s" % (parm.center(12))
     1158                ptstr += ptfmt % (Inst[parm][1])
     1159                if parm in sigDict:
     1160                    sigstr += ptfmt % (sigDict[parm])
     1161                else:
     1162                    sigstr += 12*' '
     1163        print ptlbls
     1164        print ptstr
     1165        print sigstr
     1166       
     1167    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
     1168        parmDict.update(zip(varyList,values))
     1169        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
     1170
     1171    peakPos = []
     1172    peakDsp = []
     1173    peakWt = []
     1174    for peak in IndexPeaks:
     1175        if peak[2] and peak[3]:
     1176            peakPos.append(peak[0])
     1177            peakDsp.append(peak[8])
     1178            peakWt.append(1/peak[1])
     1179    peakPos = np.array(peakPos)
     1180    peakDsp = np.array(peakDsp)
     1181    peakWt = np.array(peakWt)
     1182    dataType,insDict,insVary = SetInstParms()
     1183    parmDict = {}
     1184    parmDict.update(insDict)
     1185    varyList = insVary
     1186    while True:
     1187        begin = time.time()
     1188        values =  np.array(Dict2Values(parmDict, varyList))
     1189        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.0001,
     1190            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
     1191        ncyc = int(result[2]['nfev']/2)
     1192        runtime = time.time()-begin   
     1193        chisq = np.sum(result[2]['fvec']**2)
     1194        Values2Dict(parmDict, varyList, result[0])
     1195        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
     1196        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)
     1197        print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
     1198        print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)
     1199        try:
     1200            sig = np.sqrt(np.diag(result[1])*GOF)
     1201            if np.any(np.isnan(sig)):
     1202                print '*** Least squares aborted - some invalid esds possible ***'
     1203            break                   #refinement succeeded - finish up!
     1204        except ValueError:          #result[1] is None on singular matrix
     1205            print '**** Refinement failed - singular matrix ****'
     1206        return
     1207       
     1208    sigDict = dict(zip(varyList,sig))
     1209    GetInstParms(parmDict,Inst,varyList)
     1210    InstPrint(Inst,sigDict)
    11251211           
    11261212def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
     
    14111497        elif FitPgm == 'BFGS':
    14121498            print 'Other program here'
    1413             return
     1499            return {}
    14141500       
    14151501    sigDict = dict(zip(varyList,sig))
     
    14231509    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
    14241510    PeaksPrint(dataType,parmDict,sigDict,varyList)
     1511    return sigDict
    14251512
    14261513def calcIncident(Iparm,xdata):
  • trunk/GSASIIpwdGUI.py

    r1442 r1443  
    399399        dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
    400400        try:
    401             G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,inst2,data,oneCycle,controls,dlg)
     401            sigDict = G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,inst2,data,oneCycle,controls,dlg)
    402402        finally:
    403403            wx.EndBusyCursor()   
     
    947947                newpeaks.append(G2mth.setPeakparms(data,Inst2,peak[0],peak[2]))
    948948            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Peak List'),newpeaks)
    949                    
     949           
     950    def OnCalibrate(event):
     951        IndexPeaks = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Index Peak List'))
     952        if not len(IndexPeaks):
     953            G2frame.ErrorDialog('Can not calibrate','Index Peak List empty')
     954            return
     955#        Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))[0]
     956        G2pwd.DoCalibInst(IndexPeaks,data)
     957        UpdateInstrumentGrid(G2frame,data)
     958        XY = []
     959        for peak in IndexPeaks:
     960            if peak[2] and peak[3]:
     961                XY.append([peak[8],peak[0]])
     962        if len(XY):
     963            G2plt.PlotCalib(G2frame,data,XY,newPlot=True)
     964
    950965    def OnLoad(event):
    951966        '''Loads instrument parameters from a G2 .instprm file
     
    11761191                    dspLst.append([10,6])
    11771192                    instSizer.Add(waveVal,0,WACV)
    1178                     if ifHisto:
    1179                         refFlgElem.append([key,2])                   
    1180                         instSizer.Add(RefineBox(key),0,WACV)
    1181                     else:
    1182                         refFlgElem.append(None)                   
    1183                         instSizer.Add((5,5),0)
     1193                    refFlgElem.append([key,2])                   
     1194                    instSizer.Add(RefineBox(key),0,WACV)
     1195#                    if ifHisto:
     1196#                        refFlgElem.append([key,2])                   
     1197#                        instSizer.Add(RefineBox(key),0,WACV)
     1198#                    else:
     1199#                        refFlgElem.append(None)                   
     1200#                        instSizer.Add((5,5),0)
    11841201                for item in ['Zero','Polariz.']:
    11851202                    if item in insDef:
     
    11921209                        itemVal = G2gd.ValidatedTxtCtrl(G2frame.dataDisplay,insVal,item,nDig=(10,4),typeHint=float,OnLeave=AfterChange)
    11931210                        instSizer.Add(itemVal,0,WACV)
    1194                         if ifHisto:
    1195                             refFlgElem.append([item,2])
    1196                             instSizer.Add(RefineBox(item),0,WACV)
    1197                         else:
    1198                             refFlgElem.append(None)                   
    1199                             instSizer.Add((5,5),0)
     1211                        refFlgElem.append([item,2])
     1212                        instSizer.Add(RefineBox(item),0,WACV)
     1213#                        if ifHisto:
     1214#                            refFlgElem.append([item,2])
     1215#                            instSizer.Add(RefineBox(item),0,WACV)
     1216#                        else:
     1217#                            refFlgElem.append(None)                   
     1218#                            instSizer.Add((5,5),0)
    12001219                    else:                           #skip Polariz. for neutrons
    12011220                        instSizer.Add((5,5),0)
     
    12581277                    elemKeysLst.append([item,1])
    12591278                    dspLst.append(nDig)
    1260                     if not ifHisto and item in ['difC','difA','difB','Zero',]:
    1261                         refFlgElem.append(None)
    1262                         instSizer.Add((5,5),0)
    1263                     else:
    1264                         refFlgElem.append([item,2])
    1265                         instSizer.Add(RefineBox(item),0,WACV)
     1279                    refFlgElem.append([item,2])
     1280                    instSizer.Add(RefineBox(item),0,WACV)
     1281#                    if not ifHisto and item in ['difC','difA','difB','Zero',]:
     1282#                        refFlgElem.append(None)
     1283#                        instSizer.Add((5,5),0)
     1284#                    else:
     1285#                        refFlgElem.append([item,2])
     1286#                        instSizer.Add(RefineBox(item),0,WACV)
    12661287        elif 'S' in insVal['Type']:                       #single crystal data
    12671288            if 'C' in insVal['Type']:               #constant wavelength
     
    13571378        G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.InstMenu)
    13581379        if not G2frame.dataFrame.GetStatusBar():
    1359             Status = G2frame.dataFrame.CreateStatusBar()
     1380            Status = G2frame.dataFrame.CreateStatusBar()           
     1381        G2frame.Bind(wx.EVT_MENU,OnCalibrate,id=G2gd.wxID_INSTCALIB)
    13601382        G2frame.Bind(wx.EVT_MENU,OnLoad,id=G2gd.wxID_INSTLOAD)
    13611383        G2frame.Bind(wx.EVT_MENU,OnSave,id=G2gd.wxID_INSTSAVE)
     
    18641886            else:
    18651887                G2frame.dataDisplay.SetReadOnly(r,c,isReadOnly=True)
    1866             if 'PNT' in Inst['Type'][0] and data[r][3]:
    1867                 X = G2lat.Dsp2pos(Inst,data[r][8])
    1868                 Y = data[r][0]
    1869                 XY.append([X,Y-X])
     1888        if data[r][2] and data[r][3]:
     1889            XY.append([data[r][8],data[r][0]])
    18701890    G2frame.dataDisplay.Bind(wg.EVT_GRID_CELL_LEFT_CLICK, RefreshIndexPeaksGrid)
    18711891    G2frame.dataDisplay.Bind(wx.EVT_KEY_DOWN, KeyEditPickGrid)                 
     
    18731893    G2frame.dataDisplay.AutoSizeColumns(False)
    18741894    G2frame.dataFrame.setSizePosLeft([490,300])
    1875     if len(XY): #NB: only for TOF
    1876         G2plt.PlotXY(G2frame,np.array(XY),XY2=[],labelX='T-calc',labelY='Tobs-Tcalc',
    1877             newPlot=True,Title='Diffractometer const')
     1895    if len(XY):
     1896        G2plt.PlotCalib(G2frame,Inst,XY,newPlot=True)
    18781897     
    18791898################################################################################
Note: See TracChangeset for help on using the changeset viewer.