add pdfstream

scripts/_get_pdf.py

@@ -0,0 +1,122 @@
+from pdfstream.transformation.cli import transform
+import numpy as np
+import os
+import glob
+import matplotlib.pyplot as plt
+from scipy.signal import savgol_filter
+from pdfstream.io import load_array
+from pdfstream.transformation.io import load_pdfconfig, write_pdfgetter
+from pdfstream.transformation.main import get_pdf
+from pathlib import Path
+from diffpy.pdfgetx import PDFConfig
+
+def globfilename(path, search_item='*.tiff'):
+    filenames = glob.glob(path+search_item)
+    inputnames=[filename.replace('\\','/') for filename in filenames]
+    return(inputnames)
+
+def transform_bkg(
+    cfg_file,
+    data_file: str,
+    output_dir: str = ".",
+    plot_setting: typing.Union[str, dict] = None,
+    test: bool = False,
+) -> typing.Dict[str, str]:
+    """Transform the data."""
+    if isinstance(cfg_file,str):
+        pdfconfig = load_pdfconfig(cfg_file)
+    else:
+        pdfconfig = cfg_file
+
+    chi = load_array(data_file)
+    pdfgetter = get_pdf(pdfconfig, chi, plot_setting=plot_setting)
+    filename = Path(data_file).stem
+    dct = write_pdfgetter(output_dir, filename, pdfgetter)
+    if not test:
+        plt.show()
+    return dct
+
+
+def count_header(filename):
+    # count header 
+    with open (filename,'r') as myfile:
+        data = myfile.read().splitlines()
+    count = 0
+    for i, lines in enumerate(data):
+        if lines == '#### start data':
+            count =i
+    return count+2
+
+def find_qmax(filename, ave_cutoff=8e-03, window_length=21):
+
+    header_count=count_header(filename)
+    q, f = np.loadtxt(filename, skiprows= header_count, unpack=True)
+    qmax = q[-1]
+    fhat = savgol_filter(f, window_length, 3)
+    fdiff=f-fhat
+    zero_cross = np.where(np.diff(np.signbit(fhat)))[0]
+    z=zero_cross[::-1]
+    for i, i_zero in enumerate(z):
+        ave=np.average(abs(fdiff[i_zero:]))
+
+        if ave > ave_cutoff:
+            print(i_zero, q[i_zero], ave)
+        else:
+            qmax = np.round(q[z[i-1]],1)
+            print(i_zero, q[i_zero], ave)
+            print(f'qmax ={qmax}')
+            break    
+    plt.figure(figsize=(16,8))            
+    plt.plot(q,f)
+    plt.plot(q, fhat, 'r')
+
+    plt.plot(q, (f-fhat))   
+    print(qmax)
+    return qmax
+
+
+if __name__ == "__main__":
+
+    homepath='/nsls2/users/hzhong/pdfstream/'
+    testpath = './chi_files/'
+    output_dir=homepath
+    testfiles = globfilename(testpath, "*.chi")
+    testfile =testfiles[1]
+    bkg_file = testfiles[0]
+
+    pdfconfig = PDFConfig()
+    #pdfconfig.readConfig(cfg_file)
+    config_file= globfilename('./','*.cfg')[0]
+    pdfconfig.readConfig(config_file)
+    pdfconfig.backgroundfiles = bkg_file
+    pdfconfig.bgscale=0.9
+    print(pdfconfig)
+
+    transform('pdfconfig.cfg', testfile, output_dir=output_dir, plot_setting={'marker':'.','color':'green'} )
+    transform_bkg('pdfconfig.cfg', testfile, output_dir=output_dir, plot_setting={'marker':'.','color':'green'} )
+    transform_bkg(pdfconfig, testfile, output_dir=output_dir, plot_setting={'marker':'.','color':'green'} )
+
+    pdfconfig = PDFConfig()
+    config_file= globfilename('./','*.cfg')[0]
+    pdfconfig.readConfig(config_file)
+    pdfconfig.qmaxinst = 24.0
+    pdfconfig.qmin=1.0
+    pdfconfig.qmax=20
+    pdfconfig.rpoly=0.9
+    pdfconfig.rmin=0.1
+    pdfconfig.rmax=15
+    # pdfconfig.rpoly=0.9
+    pdfconfig.bgscale=1
+    pdfconfig.outputtypes=['sq','fq','gr']
+    print(pdfconfig)
+
+    transform_bkg(pdfconfig, testfile, output_dir=output_dir, plot_setting={'marker':'.','color':'green'} )
+
+    fqfile = globfilename(output_dir, '*/*.fq')[1]
+    print(fqfile)
+
+    qmax=find_qmax(fqfile,  ave_cutoff=8e-03, window_length=101)
+
+    pdfconfig.qmax=qmax
+    print(pdfconfig)
+    transform_bkg(pdfconfig, testfile, output_dir=output_dir, plot_setting={'marker':'.','color':'green'} )

scripts/_pdf_calculator.py

@@ -20,7 +20,7 @@ def _no_oxidation_cif(cif_file):
 
 
 
-def _pdffit2_CsPbX3(gr_data, cif_list, qmax=18, qdamp=0.031, qbroad=0.032, fix_APD=True, toler=0.000001):
+def _pdffit2_CsPbX3(gr_data, cif_list, qmax=18, qdamp=0.031, qbroad=0.032, fix_APD=True, toler=0.000001, return_pf=False):
 
     # Initialize the CifParser with the path to your .cif file
     # Parse the .cif file
@@ -57,7 +57,7 @@ def _pdffit2_CsPbX3(gr_data, cif_list, qmax=18, qdamp=0.031, qbroad=0.032, fix_A
     pf.setvar(pf.qbroad, qbroad)
 
     # Refine 
-    pf.pdfrange(1, 2.5, 60)
+    pf.pdfrange(1, 2.5, 80)
     pf.refine(toler=toler)
 
 
@@ -68,7 +68,10 @@ def _pdffit2_CsPbX3(gr_data, cif_list, qmax=18, qdamp=0.031, qbroad=0.032, fix_A
         pf.setphase(i+1)
         particel_size.append(pf.getvar(pf.spdiameter))
 
-    return phase_fraction, particel_size
+    if return_pf:
+        return pf
+    else:
+        return phase_fraction, particel_size
 
 
 
@@ -90,16 +93,18 @@ def _set_CsPbBr3_constrain(PDF_calculator_object, phase_idx=1, fix_APD=True):
 
     # Refine phase scale factor.  Right side can have formulas.
     pf.constrain('pscale', '@111')
-    pf.setpar(111, 0.5)
+    pf.setpar(111, 0.85)
     # pf.setpar(20, pf.getvar(pf.pscale) / 2.0)
 
     # Refine sharpening factor for correlated motion of close atoms.
     pf.constrain(pf.delta2, '@122')
     pf.setpar(122, 6.87)
+    pf.fixpar(122)
 
     # Refine diameter for the spherical particle
     pf.constrain(pf.spdiameter, '@133')
-    pf.setpar(133, 46)
+    pf.setpar(133, 80)
+    # pf.fixpar(133)
 
     # Set temperature factors isotropic to each atom
     # idx starts from 1 not 0