Changeset 4519 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Jul 13, 2020 1:27:06 PM (15 months ago)
Author:
vondreele
Message:

G2fpaGUI - minor cleanup
G2lattice - test on 'T' instead of 'C' so that 'B' PWDR types properly handled
G2math - add getPinkalpha, getPinkbeta % deriv routines for pink beam peak shapes

modify setPeakparams to use them

G2plot - test on 'T' instead of 'C' so that 'B' PWDR types properly handled
G2pwd - add pink beam peak shape function - same as TOF peak shape; scale sig & gam to be in centidegrees

  • add wtFactor to FitPeaks? fix bad PWDR weights for peak fitting so reduced chi2 is more appropriate

G2pwdGUI - test on 'T' instead of 'C' so that 'B' PWDR types properly handled

  • add wtFactor to DoPeakFit? to fix bad PWDR weights for peak fitting so reduced chi2 is more appropriate

G2strMath - latest in incommensurate mag str factors - no real improvement

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r4518 r4519  
    801801    alpTOF = lambda dsp,alp: alp/dsp
    802802    betTOF = lambda dsp,bet0,bet1,betq: bet0+bet1/dsp**4+betq/dsp**2
    803     if 'C' in Inst['Type'][0]:
    804         s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
    805         g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
    806         return getgamFW(g,s)/100.  #returns FWHM in deg
    807     else:
     803    if 'T' in Inst['Type'][0]:
    808804        dsp = pos/Inst['difC'][0]
    809805        alp = alpTOF(dsp,Inst['alpha'][0])
     
    812808        g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
    813809        return getgamFW(g,s)+np.log(2.0)*(alp+bet)/(alp*bet)
     810    else:
     811        s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
     812        g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
     813        return getgamFW(g,s)/100.  #returns FWHM in deg
    814814   
    815815def getgamFW(g,s):
     
    854854    if 'T' in dataType:
    855855        q = 2.*np.pi*parmDict[pfx+'difC']/xdata
    856     elif 'C' in dataType:
     856    else:
    857857        wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
    858858        q = npT2q(xdata,wave)
     
    990990    if 'T' in dataType:
    991991        q = 2.*np.pi*parmDict[hfx+'difC']/xdata
    992     elif 'C' in dataType:
     992    else:
    993993        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
    994994        q = 2.*np.pi*npsind(xdata/2.)/wave
     
    10871087            if 'C' in dataType:
    10881088                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
    1089             else: #'T'OF
     1089            else: #'T' or 'B'
    10901090                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
    10911091            iBeg = np.searchsorted(xdata,pkP-fmin)
     
    13211321            except KeyError:        #no more peaks to process
    13221322                return yb+yc
     1323    elif 'B' in dataType:
     1324        iPeak = 0
     1325        dsp = 1.0 #for now - fix later
     1326        while True:
     1327            try:
     1328                pos = parmDict['pos'+str(iPeak)]
     1329                tth = (pos-parmDict['Zero'])
     1330                intens = parmDict['int'+str(iPeak)]
     1331                alpName = 'alp'+str(iPeak)
     1332                if alpName in varyList:
     1333                    alp = parmDict[alpName]
     1334                else:
     1335                    alp = G2mth.getPinkalpha(parmDict,tth)
     1336                alp = max(0.1,alp)
     1337                betName = 'bet'+str(iPeak)
     1338                if betName in varyList:
     1339                    bet = parmDict[betName]
     1340                else:
     1341                    bet = G2mth.getPinkbeta(parmDict,tth)
     1342                bet = max(0.0001,bet)
     1343                sigName = 'sig'+str(iPeak)
     1344                if sigName in varyList:
     1345                    sig = parmDict[sigName]
     1346                else:
     1347                    sig = G2mth.getCWsig(parmDict,tth)
     1348                sig = max(sig,0.001)          #avoid neg sigma^2
     1349                gamName = 'gam'+str(iPeak)
     1350                if gamName in varyList:
     1351                    gam = parmDict[gamName]
     1352                else:
     1353                    gam = G2mth.getCWgam(parmDict,tth)
     1354                gam = max(gam,0.001)             #avoid neg gamma
     1355                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
     1356                iBeg = np.searchsorted(xdata,pos-fmin)
     1357                iFin = np.searchsorted(xdata,pos+fmin)
     1358                if not iBeg+iFin:       #peak below low limit
     1359                    iPeak += 1
     1360                    continue
     1361                elif not iBeg-iFin:     #peak above high limit
     1362                    return yb+yc
     1363                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])
     1364                iPeak += 1
     1365            except KeyError:        #no more peaks to process
     1366                return yb+yc       
    13231367    else:
    13241368        Pdabc = parmDict['Pdabc']
     
    14801524            except KeyError:        #no more peaks to process
    14811525                break
     1526    elif 'B' in dataType:
     1527        iPeak = 0
     1528        while True:
     1529            try:
     1530                pos = parmDict['pos'+str(iPeak)]
     1531                tth = (pos-parmDict['Zero'])
     1532                intens = parmDict['int'+str(iPeak)]
     1533                alpName = 'alp'+str(iPeak)
     1534                if alpName in varyList:
     1535                    alp = parmDict[alpName]
     1536                else:
     1537                    alp = G2mth.getPinkalpha(parmDict,tth)
     1538                    dada0,dada1 = G2mth.getPinkalphaDeriv(tth)
     1539                alp = max(0.0001,alp)
     1540                betName = 'bet'+str(iPeak)
     1541                if betName in varyList:
     1542                    bet = parmDict[betName]
     1543                else:
     1544                    bet = G2mth.getPinkbeta(parmDict,tth)
     1545                    dbdb0,dbdb1 = G2mth.getPinkbetaDeriv(tth)
     1546                bet = max(0.0001,bet)
     1547                sigName = 'sig'+str(iPeak)
     1548                if sigName in varyList:
     1549                    sig = parmDict[sigName]
     1550                    dsdU = dsdV = dsdW = 0
     1551                else:
     1552                    sig = G2mth.getCWsig(parmDict,tth)
     1553                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
     1554                sig = max(sig,0.001)          #avoid neg sigma
     1555                gamName = 'gam'+str(iPeak)
     1556                if gamName in varyList:
     1557                    gam = parmDict[gamName]
     1558                    dgdX = dgdY = dgdZ = 0
     1559                else:
     1560                    gam = G2mth.getCWgam(parmDict,tth)
     1561                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
     1562                gam = max(gam,0.001)             #avoid neg gamma
     1563                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig/1.e4,gam/100.)
     1564                iBeg = np.searchsorted(xdata,pos-fmin)
     1565                iFin = np.searchsorted(xdata,pos+fmin)
     1566                if not iBeg+iFin:       #peak below low limit
     1567                    iPeak += 1
     1568                    continue
     1569                elif not iBeg-iFin:     #peak above high limit
     1570                    break
     1571                dMdpk = np.zeros(shape=(7,len(xdata)))
     1572                dMdipk = getdEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])
     1573                for i in range(1,6):
     1574                    dMdpk[i][iBeg:iFin] += cw[iBeg:iFin]*intens*dMdipk[i]
     1575                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
     1576                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4]/1.e4,'gam':dMdpk[5]/100.}
     1577                for parmName in ['pos','int','alp','bet','sig','gam']:
     1578                    try:
     1579                        idx = varyList.index(parmName+str(iPeak))
     1580                        dMdv[idx] = dervDict[parmName]
     1581                    except ValueError:
     1582                        pass
     1583                if 'U' in varyList:
     1584                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
     1585                if 'V' in varyList:
     1586                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
     1587                if 'W' in varyList:
     1588                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
     1589                if 'X' in varyList:
     1590                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
     1591                if 'Y' in varyList:
     1592                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
     1593                if 'Z' in varyList:
     1594                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
     1595                if 'alpha-0' in varyList:
     1596                    dMdv[varyList.index('alpha-0')] += dada0*dervDict['alp']
     1597                if 'alpha-1' in varyList:
     1598                    dMdv[varyList.index('alpha-1')] += dada1*dervDict['alp']
     1599                if 'beta-0' in varyList:
     1600                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
     1601                if 'beta-1' in varyList:
     1602                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
     1603                iPeak += 1
     1604            except KeyError:        #no more peaks to process
     1605                break       
    14821606    else:
    14831607        Pdabc = parmDict['Pdabc']
     
    17201844    return True
    17211845           
    1722 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
     1846def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,wtFactor=1.0,dlg=None):
    17231847    '''Called to perform a peak fit, refining the selected items in the peak
    17241848    table as well as selected items in the background.
     
    17441868    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
    17451869      and derivType = controls['deriv type']. If None default values are used.
     1870    :param float wtFactor. weight multiplier; = 1.0 by default
    17461871    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
    17471872      Defaults to None, which means no updates are done.
     
    18421967            insVals.append(Inst[parm][1])
    18431968            if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
    1844                 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
     1969                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1'] and Inst[parm][2]:
    18451970                    insVary.append(parm)
    18461971        instDict = dict(zip(insNames,insVals))
     
    18601985                pos = parmDict['pos'+str(iPeak)]
    18611986                if sigName not in varyList:
    1862                     if 'C' in Inst['Type'][0]:
    1863                         parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
    1864                     else:
     1987                    if 'T' in Inst['Type'][0]:
    18651988                        dsp = G2lat.Pos2dsp(Inst,pos)
    18661989                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
     1990                    else:
     1991                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
    18671992                gamName = 'gam'+str(iPeak)
    18681993                if gamName not in varyList:
    1869                     if 'C' in Inst['Type'][0]:
    1870                         parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
    1871                     else:
     1994                    if 'T' in Inst['Type'][0]:
    18721995                        dsp = G2lat.Pos2dsp(Inst,pos)
    18731996                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
     1997                    else:
     1998                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
    18741999                iPeak += 1
    18752000            except KeyError:
     
    18842009        for parm in Inst:
    18852010            if parm in  ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
    1886                 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
     2011                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1']:
    18872012                ptlbls += "%s" % (parm.center(12))
    18882013                ptstr += ptfmt % (Inst[parm][1])
     
    19012026        if 'C' in dataType:
    19022027            names = ['pos','int','sig','gam']
    1903         else:
     2028        else:   #'T' and 'B'
    19042029            names = ['pos','int','alp','bet','sig','gam']
    19052030        for i,peak in enumerate(Peaks):
     
    19152040        if 'C' in Inst['Type'][0]:
    19162041            names = ['pos','int','sig','gam']
    1917         else:   #'T'
     2042        else:   #'T' & 'B'
    19182043            names = ['pos','int','alp','bet','sig','gam']
    19192044        for i,peak in enumerate(Peaks):
     
    19252050                if parName in varyList:
    19262051                    peak[2*j] = parmDict[parName]
    1927                 elif 'alpha' in parName:
    1928                     peak[2*j] = parmDict['alpha']/dsp
    1929                 elif 'beta' in parName:
    1930                     peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
     2052                elif 'alp' in parName:
     2053                    if 'T' in Inst['Type'][0]:
     2054                        peak[2*j] = G2mth.getTOFalpha(parmDict,dsp)
     2055                    else: #'B'
     2056                        peak[2*j] = G2mth.getPinkalpha(parmDict,pos)
     2057                elif 'bet' in parName:
     2058                    if 'T' in Inst['Type'][0]:
     2059                        peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
     2060                    else:   #'B'
     2061                        peak[2*j] = G2mth.getPinkbeta(parmDict,pos)
    19312062                elif 'sig' in parName:
    1932                     if 'C' in Inst['Type'][0]:
     2063                    if 'T' in Inst['Type'][0]:
     2064                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
     2065                    else:   #'C' & 'B'
    19332066                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
    1934                     else:
    1935                         peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
    19362067                elif 'gam' in parName:
    1937                     if 'C' in Inst['Type'][0]:
     2068                    if 'T' in Inst['Type'][0]:
     2069                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
     2070                    else:   #'C' & 'B'
    19382071                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
    1939                     else:
    1940                         peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
    19412072                       
    19422073    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
     
    19442075        if 'C' in dataType:
    19452076            names = ['pos','int','sig','gam']
    1946         else:   #'T'
     2077        else:   #'T' & 'B'
    19472078            names = ['pos','int','alp','bet','sig','gam']           
    19482079        head = 13*' '
     
    19562087        if 'C' in dataType:
    19572088            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
    1958         else:
     2089        else:   #'T' & 'B'
    19592090            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
    19602091        for i,peak in enumerate(Peaks):
     
    20242155        badVary = []
    20252156        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
    2026                args=(x[xBeg:xFin],(y+fixback)[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
     2157               args=(x[xBeg:xFin],(y+fixback)[xBeg:xFin],wtFactor*w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
    20272158        ncyc = int(result[2]['nfev']/2)
    20282159        runtime = time.time()-begin   
    20292160        chisq = np.sum(result[2]['fvec']**2)
    20302161        Values2Dict(parmDict, varyList, result[0])
    2031         Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*(y+fixback)[xBeg:xFin]**2))*100.      #to %
     2162        Rvals['Rwp'] = np.sqrt(chisq/np.sum(wtFactor*w[xBeg:xFin]*(y+fixback)[xBeg:xFin]**2))*100.      #to %
    20322163        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
    20332164        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))
Note: See TracChangeset for help on using the changeset viewer.