notebook.community

Edit and run

XPCS Pipeline

"This notebook corresponds to version {{ version }} of the pipeline tool: https://github.com/NSLS-II/pipelines"

NSLS2 data retrieval imports





In [1]:



    


from databroker import DataBroker as db, get_images, get_table, get_events
from filestore.api import register_handler, deregister_handler
from filestore.retrieve import _h_registry, _HANDLER_CACHE



    











In [ ]:



    


hdr = db[{{ uid }}]



    











In [2]:



    


import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import time
from ipywidgets import interact



    











In [3]:



    


cd /home/yuzhang/chx-pipelines/Develops/



    


















    






/home/yuzhang/chx-pipelines/Develops

















In [ ]:



    


cd /XF11ID/analysis/Analysis_Pipelines/Develop/



    











In [4]:



    


%run develop.py
%run two_time.py



    


















    






/home/yuzhang/.conda/envs/user_analysis/lib/python3.4/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.
  "`IPython.html.widgets` has moved to `ipywidgets`.", ShimWarning)

















In [5]:



    


%matplotlib notebook
#%matplotlib inline

Users put

  • uid for Bluesky Scan
  • filename for acquiring data directly by EigerSofteare




In [6]:



    


BlueScan = True
DirectAcq = False
detector = 'eiger_4M_cam_img_image_lightfield'  #for 4M

Users put uid here





In [22]:



    


if BlueScan:
    uid = '54614d43'
    #uid = '95782687'
    uid = '95782687'
    uid= 'ff9f20c0'
    uid='71720966'
    uid='1663d34a'
    uid = 'f505e052-3baa-47d4-bdc4-61c2eb1bcc7a'  #sid= 551, 1%PEG, 
    uid='ee6975a1-9161'   #1% wt PEG
    uid='ec41569e'  #a gisaxs example, 
    uid = 'ffe9d518'  # 10 mTorr 1sec/frame
 
    
else:
    uid = '/XF11ID/data/2015/11/23/d01ab510-3cf3-4719-bee3_795_master.h5'

Get data





In [23]:



    


if BlueScan:
    hdr = db[uid]
    ev, = get_events(  hdr, [detector] )
    imgs = ev['data'][detector]
else:    
    imgs =  Images(uid)
print (imgs)
Nimg=len(imgs)



    


















    






hdf5 path = /XF11ID/data/2015/11/8/3f84f57d-79b7-4a10-8b51_5_master.h5
<Frames>
Length: 15000 frames
Frame Shape: 2167 x 2070
Pixel Datatype: uint32

Get data path





In [24]:



    


if BlueScan:
    from datetime import datetime
    dt = datetime.fromtimestamp(hdr['start'].time)
    path ='/XF11ID/analysis' + '/%s/%s/%s/' % (dt.year, dt.month, dt.day)
else:
    path ='/XF11ID/analysis/2015/11/23/' 
path



    


















    
Out[24]:








'/XF11ID/analysis/2015/11/8/'

Note: experiment information

  • The physical size of the pixels
  • Wavelegth of the X-rays - (units in Angstroms)
  • Detector to sample distance
  • Exposure time - (units in seconds)
  • acqusition period - (units in seconds)
  • dead time - (units in seconds)
  • time per frame = (exposure time + dead_time or acqusition period) - (units in seconds)




In [25]:



    


imgs.md



    


















    
Out[25]:








{'beam_center_x': 840.0,
 'beam_center_y': 1830.0,
 'count_time': 1.0,
 'detector_distance': 4.8400002,
 'frame_time': 1.00001,
 'framerate': 0.99998998651994853,
 'incident_wavelength': 1.3776,
 'pixel_mask': array([[1, 1, 1, ..., 1, 1, 0],
        [1, 1, 1, ..., 1, 1, 1],
        [1, 1, 1, ..., 1, 1, 1],
        ..., 
        [1, 1, 1, ..., 1, 1, 1],
        [1, 1, 1, ..., 1, 1, 1],
        [1, 1, 1, ..., 1, 1, 1]], dtype=uint32),
 'x_pixel_size': 7.5000004e-05,
 'y_pixel_size': 7.5000004e-05}

















In [26]:



    


# The physical size of the pixels
dpix = imgs.md['x_pixel_size'] * 1000.  
lambda_ = imgs.md['incident_wavelength']    # wavelegth of the X-rays in Angstroms
Ldet = 4812.        # detector to sample distance (mm)

exposuretime= imgs.md['count_time']
acquisition_period = imgs.md['frame_time']

# deadtime= 0   # 60e-6 
# timeperframe = exposuretime + deadtime
timeperframe = acquisition_period  

timeperframe, exposuretime



    


















    
Out[26]:








(1.00001, 1.0)

load a mask

load the mask saved in the mask_pipeline





In [27]:



    


mask = np.load( path +  str(uid)+ "_mask.npy")

Reverse the mask in y-direction due to the coordination difference between python and Eiger software





In [28]:



    


maskr = mask[::-1,:]

Plot the mask





In [29]:



    


fig, ax = plt.subplots()
im=ax.imshow(maskr, origin='lower' ,vmin=0,vmax=1,cmap='viridis')
fig.colorbar(im)
plt.show()

Interactive way to browse through images.

Note : Provide the number of images that you want to browse





In [30]:



    


def view_image(i):    
    fig, ax = plt.subplots()
    ax.imshow(imgs[i]*mask, interpolation='nearest', cmap='viridis',
                  origin='lower', norm= LogNorm(vmin=0.001, vmax=1e1) )
    ax.set_title("Browse the Image Stack")
    plt.show()



    











In [31]:



    


#interact(view_image, i=(0, Nimg-1))

Movie





In [32]:



    


def view_image(sleeps=1, ims=0, ime = 1):    
    fig, ax = plt.subplots()  
    for i in range( ims, ime  ):
        im=ax.imshow(imgs[i]*mask,  interpolation='nearest', cmap='viridis',
                  origin='lower', norm= LogNorm( vmin=0.01, vmax=10 ) )
        ax.set_title("images_%s"%i)
        
        time.sleep( sleeps )
        plt.draw()
    #fig.colorbar(im)
        
#view_image(.2, 0, 2)

hey, let's see if any images are bad!

load the image intensity (kymograph) saved in the mask_pipeline





In [33]:



    


kymo_sum = np.load( path +  str(uid)+"_kymo_sum.npy" )



    











In [34]:



    


fig, axes = plt.subplots(  )
axes.plot( kymo_sum, '-go'  ) 
ax.set_ylabel('Intensity')
ax.set_xlabel('Frame')
ax.set_title('Kymograph_sum') 
plt.show()

Get the Averaged Image Data

load the average intensity saved in the mask_pipeline





In [35]:



    


avg_img = np.load( path + str(uid)+"_avg_img.npy" )
avg_imgm =  avg_img * mask

Reverse the image in y-direction due to the coordination difference between python and Eiger software





In [36]:



    


avg_imgr  = avg_img[::-1,:] 
avg_imgmr  = avg_imgm[::-1,:]

NOTE:

  • mask: the mask file in the python conventional coordination
  • maskr: the mask file in the eiger coordination (maksr = mask[::-1,:])
  • avg_img: the average intensity frames in the python conventional coordination
  • avg_imgm: the maksed average intensity frames in the python conventional coordination (avg_imgm = avg_img * mask)
  • avg_imgr: the average intensity frames in the eiger coordination (avg_imgr = avg_img[::-1,:] )
  • avg_imgmr: the maksed average intensity frames in the eiger coordination (avg_imgmr = avg_imgm[::-1,:] )

Plot the averged image with the mask





In [37]:



    


fig, ax = plt.subplots()
im = ax.imshow(avg_imgmr, cmap='viridis',origin='lower',
               norm= LogNorm(vmin=0.001, vmax=1e1))
ax.set_title("Masked Averaged Image_Reversed")
fig.colorbar(im)
plt.show()

Note: experiment information

  • The physical size of the pixels
  • Wavelegth of the X-rays - (units in Angstroms)
  • Detector to sample distance
  • Exposure time - (units in seconds)
  • acqusition period - (units in seconds)
  • dead time - (units in seconds)
  • time per frame = (exposure time + dead_time or acqusition period) - (units in seconds)




In [38]:



    


# The physical size of the pixels
dpix = imgs.md['x_pixel_size'] * 1000.  
lambda_ = imgs.md['incident_wavelength']    # wavelegth of the X-rays in Angstroms
Ldet = 4810        # detector to sample distance (mm)

exposuretime= imgs.md['count_time']
acquisition_period = imgs.md['frame_time']

# deadtime= 0   # 60e-6 
# timeperframe = exposuretime + deadtime
timeperframe = acquisition_period  

timeperframe, exposuretime



    


















    
Out[38]:








(1.00001, 1.0)

GiSAXS get q-map (qz and qr)

Users put incident-Beam and Reflection_Beam Centers here!!!





In [39]:



    


inc_x0 = 1871
inc_y0 = 339
refl_x0 = 1871
refl_y0 = 811      #1670
Lsd= 4.81
lamda= lambda_  #12.4/9



    











In [40]:



    


lamda



    


















    
Out[40]:








1.3776

















In [41]:



    


alphaf,thetaf, alphai, phi = get_reflected_angles( inc_x0, inc_y0,refl_x0 , refl_y0, Lsd=Lsd )



    


















    






The incident_angle (alphai) is: 0.210835132612

















In [42]:



    


qx, qy, qr, qz = convert_gisaxs_pixel_to_q( inc_x0, inc_y0,refl_x0,refl_y0, lamda=lamda, Lsd=Lsd )



    


















    






The incident_angle (alphai) is: 0.210835132612

















In [43]:



    


fig, ax = plt.subplots()
#im=ax.imshow(alphaf, origin='lower' ,cmap='viridis',norm= LogNorm(vmin=0.0001,vmax=2.00))
im=ax.imshow(alphaf*180/np.pi, origin='lower' ,cmap='viridis',vmin=-1,vmax= 1.5 )
fig.colorbar(im)
ax.set_title( 'alphaf')
plt.show()



    


















    




















    























In [44]:



    


fig, ax = plt.subplots()
#im=ax.imshow(alphaf, origin='lower' ,cmap='viridis',norm= LogNorm(vmin=0.0001,vmax=2.00))
im=ax.imshow(qr, origin='lower' ,cmap='viridis',vmin=qr.min(),vmax= qr.max() )
fig.colorbar(im)
ax.set_title( 'Q-Parallel')
plt.show()



    


















    




















    























In [45]:



    


fig, ax = plt.subplots()
#im=ax.imshow(alphaf, origin='lower' ,cmap='viridis',norm= LogNorm(vmin=0.0001,vmax=2.00))
im=ax.imshow(qz, origin='lower' ,cmap='viridis',vmin=qz.min(),vmax= qz.max() )
fig.colorbar(im)
ax.set_title( 'Q-z')
plt.show()

Get q||-intensity

Users put cuts here!!!





In [49]:



    


fig, ax = plt.subplots()
im = ax.imshow(avg_imgr, cmap='viridis', origin = 'lower',  norm= LogNorm( vmin=0.00001, vmax=.5e2 ) )
ax.set_title("Masked Averaged Image")
fig.colorbar(im)
plt.show()



    


















    




















    























In [70]:



    


vert_rect = ( ( 850, 0, 980- 850, 1600-0) , ( 570, 0, 700- 570, 1600-0) )    #(y,x, hight, wdith)



    











In [71]:



    


new_mask = np.ones_like( avg_imgr)
new_mask[ :, 1020:1045] =0



    











In [72]:



    


get_qr_intensity( qr, avg_imgr, vert_rect, mask=new_mask, show_roi=True)



    


















    




















    
















    




















    
















    






0 (850, 0, 130, 1600)
1 (570, 0, 130, 1600)

User gives the interested Qz,Qr here!!!





In [73]:



    


qz_start = qz[670,0]
qz_end = qz[950,0]
qz_num= 2


qr_start =  qr[600,1700]
qr_end = qr[600,0]
qr_num = 15

Create label array (Qz, Qr, Qzr boxes)





In [74]:



    


qr_edge, qr_center = get_qedge(qr_start , qr_end, ( qr_end- qr_start)/qr_num, qr_num )
qz_edge, qz_center = get_qedge( qz_start,   qz_end,   (qz_end - qz_start)/(qz_num -0) , qz_num )

label_array_qz = get_qmap_label( qz, qz_edge)
label_array_qr = get_qmap_label( qr, qr_edge)
label_array_qzr,qzc,qrc = get_qzrmap(label_array_qz, label_array_qr,qz_center, qr_center  )

labels_qzr, indices_qzr = roi.extract_label_indices( label_array_qzr  )
labels_qz, indices_qz = roi.extract_label_indices( label_array_qz  )
labels_qr, indices_qr = roi.extract_label_indices( label_array_qr  )

num_qz = len(np.unique( labels_qz ))
num_qr = len(np.unique( labels_qr ))
num_qzr = len(np.unique( labels_qzr ))



    











In [75]:



    


num_qz,num_qr,num_qzr



    


















    
Out[75]:








(2, 15, 30)

















In [76]:



    


fig, ax = plt.subplots()
#im=ax.imshow(alphaf, origin='lower' ,cmap='viridis',norm= LogNorm(vmin=0.0001,vmax=2.00))
im=ax.imshow(label_array_qz, origin='lower' ,cmap='viridis',vmin=0,vmax= None )
fig.colorbar(im)
ax.set_title( 'Q-z_label')
plt.show()



    


















    




















    























In [77]:



    


fig, ax = plt.subplots()
#im=ax.imshow(alphaf, origin='lower' ,cmap='viridis',norm= LogNorm(vmin=0.0001,vmax=2.00))
im=ax.imshow(label_array_qr, origin='lower' ,cmap='viridis',vmin=0,vmax= None )
fig.colorbar(im)
ax.set_title( 'Q-r_label')
plt.show()

Extract the labeled array





In [91]:



    


boxes = label_array_qzr 
box_maskr = boxes*maskr



    











In [92]:



    


qind, pixelist = roi.extract_label_indices(   box_maskr  )
noqs = len( np.unique(qind) )
nopr = np.bincount(qind, minlength=(noqs+1))[1:]

Number of pixels in each q ring





In [90]:



    


nopr



    


















    
Out[90]:








array([27383, 15846, 15707, 15707, 15846, 15290, 15845, 15697, 15706,
       15846, 15703, 15843, 15707, 15846,  7922, 27776, 16074, 15933,
       15933, 16074, 15508, 16074, 15933, 15933, 16039, 15933, 16074,
       15933, 16074,  8037])

plot average image with interested Q-regions (boxes)





In [87]:



    


# plot the figure
fig, axes = plt.subplots(figsize=(8,8))
axes.set_title("Labeled Array on Averaged Data")
im,im_label = show_label_array_on_image(axes, avg_imgr, box_maskr, imshow_cmap='viridis',
                        cmap='Paired',
                         vmin=0.01, vmax=30. ,  origin="lower")
#rwidth = 200 
#x1,x2 = [center[1] - rwidth, center[1] + rwidth]
#y1,y2 = [center[0] - rwidth, center[0] + rwidth]
#axes.set_xlim( [x1,x2])
#axes.set_ylim( [y1,y2])

#fig.colorbar(im)
fig.colorbar(im_label)
plt.show()

Kymograph(waterfall plot) of the max-intensity ring





In [88]:



    


imgs_ =imgs
imgsr = Reverse_Coordinate(imgs_, maskr)



    











In [93]:



    


t0 = time.time()
data_pixel =   Get_Pixel_Array( imgsr, pixelist).get_data()
run_time(t0)



    


















    






Total time: 4.70 min

users put the number of ring with max intensity here





In [113]:



    


max_inten_ring =2



    











In [114]:



    


#kymo = roi.kymograph(imgsr, ring_mask, num = max_inten_ring)



    











In [115]:



    


pixelist_qi =  np.where( qind == max_inten_ring )[0]         
data_pixel_qi = data_pixel[:,pixelist_qi]



    











In [121]:



    


data_pixel_qi.shape



    


















    
Out[121]:








(15000, 15846)

















In [120]:



    


fig, ax = plt.subplots(figsize=(8,6))
ax.set_ylabel('Pixel')
ax.set_xlabel('Frame')
ax.set_title('Kymograph')

im = ax.imshow(data_pixel_qi.T, cmap='viridis', vmax=1.0)
fig.colorbar( im   )
ax.set_aspect(1.)
plt.show()



    


















    




















    























In [132]:



    


fig, ax = plt.subplots(figsize=(8,6))
ax.set_ylabel('Pixel')
ax.set_xlabel('Frame')
ax.set_title('Kymograph_4000-6000')

im = ax.imshow(data_pixel_qi[4000:6000,:].T, cmap='viridis', vmax=1.0)
fig.colorbar( im   )
ax.set_aspect(.1)
plt.show()

Mean intensities for each ring





In [101]:



    


mean_inten = get_mean_intensity( data_pixel, qind)



    











In [200]:



    


times = np.arange(  mean_inten[1].shape[0]   )  #*timeperframe  # get the time for each frame

fig, ax = plt.subplots(figsize=(8, 6))
ax.set_title("Mean intensity of each ring")
for i in range(  num_qzr  ):
    ax.plot(times, mean_inten[i+1], '--o', label="Ring "+str(i+1))
    ax.set_xlabel("Frame")
    ax.set_ylabel("Mean Intensity")
    ax.set_ylim(0, 6e4)
    
ax.legend( fontsize = 6) 
plt.show()

One time Correlation

Note : Enter the number of levels and number of buffers for Muliti tau one time correlation number of buffers has to be even. More details in https://github.com/scikit-xray/scikit-xray/blob/master/skxray/core/correlation.py





In [133]:



    


nopr



    


















    
Out[133]:








array([27383, 15846, 15707, 15707, 15846, 15290, 15845, 15697, 15706,
       15846, 15703, 15843, 15707, 15846,  7922, 27776, 16074, 15933,
       15933, 16074, 15508, 16074, 15933, 15933, 16039, 15933, 16074,
       15933, 16074,  8037])

















In [226]:



    


%run two_time.py



    











In [224]:



    


good_start = 8000
good_end = 9201 

imgs_ =imgs[good_start: good_end-1]
imgsr = Reverse_Coordinate(imgs_, mask) 

num_buf = 8
g2_tx1, lag_steps_tx1 = autocor_one_time( num_buf,  box_maskr, imgsr, num_lev=None,
                                bad_images=None, threshold= None )



    


















    






The lev number is 9
Doing g2 caculation of 1200 frames---
##########Total time: 3.08 min

















In [229]:



    


good_start = 8000
good_end = 9201 

imgs_ =imgs[good_start: good_end-1]
imgsr = Reverse_Coordinate(imgs_, mask) 

num_buf = 8
g2_tx2, lag_steps_tx2 = autocor_one_time( num_buf,  box_maskr, imgsr, num_lev=None,
                                bad_images=[10,150,470,570,890,], threshold= None )



    


















    






The lev number is 9
Doing g2 caculation of 1200 frames---
#bad image: 10 here!
#bad image: 150 here!
##bad image: 470 here!
#bad image: 570 here!
###bad image: 890 here!
##{1: [5, 10, 10, 10, 10, 10, 10, 10], 2: [10, 10, 9, 9], 3: [9, 9, 9, 9], 4: [9, 9, 9, 9], 5: [9, 9, 8, 9], 6: [8, 8, 8, 8], 7: [8, 5, 5, 6], 8: [5, 3, 2, 2], 9: [0, 0, 0, 0]}
Total time: 3.00 min

















In [231]:



    


good_start = 8000
good_end = 9201 

imgs_ =imgs[good_start: good_end-1]
imgsr = Reverse_Coordinate(imgs_, mask) 

num_buf = 8
g2_tx3, lag_steps_tx3 = autocor_one_time( num_buf,  box_maskr, imgsr, num_lev=None,
                                bad_images=[570], threshold= None )



    


















    






The lev number is 9
Doing g2 caculation of 1200 frames---
#####bad image: 570 here!
#####{1: [1, 2, 2, 2, 2, 2, 2, 2], 2: [2, 2, 2, 2], 3: [2, 2, 2, 2], 4: [2, 2, 2, 2], 5: [2, 2, 2, 2], 6: [2, 2, 2, 2], 7: [2, 2, 2, 2], 8: [2, 0, 0, 0], 9: [0, 0, 0, 0]}
Total time: 3.01 min

















In [238]:



    


fig,ax = plt.subplots()
ax.semilogx( lag_steps_tx3*timeperframe, g2_tx3[:,1], '--bs', label='inv-1_bad_images')
ax.semilogx( lag_steps_tx2*timeperframe, g2_tx2[:,1], '--ko', label='inv-5_bad_images')
ax.semilogx( lag_steps_tx1*timeperframe, g2_tx1[:,1], '-r', label='all_good_images')
ax.set_ylim([ 1.12,1.16 ])
ax.legend()
ax.set_title('test_dealing_with_bad_frames_in_multi_tau_code', fontsize=14,y =1.03) 
ax.set_xlabel(r"$\tau $ $(s)$", fontsize=16) 
ax.set_ylabel(r"$g_2$", fontsize=16)



    


















    




















    
















    
Out[238]:








<matplotlib.text.Text at 0x7f9fe8ebe6d8>

Another check for algorithm dealing with bad-frames





In [259]:



    


%run two_time.py



    











In [267]:



    


good_start = 7245
good_end = 9246

imgs_ =imgs[good_start: good_end-1]
imgsr = Reverse_Coordinate(imgs_, mask) 

num_buf = 8
g2_ty2, lag_steps_ty2 = autocor_one_time( num_buf,  box_maskr, imgsr, num_lev=None,
                                bad_images=None, threshold= 5000 )



    


















    






The lev number is 10
Doing g2 caculation of 2000 frames---
bad image: 0 here!
##########Total time: 4.96 min

















In [268]:



    


good_start = 7246
good_end = 9246

imgs_ =imgs[good_start: good_end-1]
imgsr = Reverse_Coordinate(imgs_, mask) 

num_buf = 8
g2_ty1, lag_steps_ty1= autocor_one_time( num_buf,  box_maskr, imgsr, num_lev=None,
                                bad_images=None, threshold= None)



    


















    






The lev number is 10
Doing g2 caculation of 1999 frames---
###########Total time: 4.95 min

















In [271]:



    


fig,ax = plt.subplots()
ax.semilogx( lag_steps_ty2*timeperframe, g2_ty2[:,1], '--bs', label='inv-1_bad_images')
ax.semilogx( lag_steps_ty1*timeperframe, g2_ty1[:,1], '--ko', label='all_good_images')
#ax.semilogx( lag_steps_tx1*timeperframe, g2_tx1[:,1], '-r', label='all_good_images')
ax.set_ylim([ 1.10,1.16 ])
ax.legend()
ax.set_title('test_dealing_with_bad_frames_in_multi_tau_code', fontsize=14,y =1.03) 
ax.set_xlabel(r"$\tau $ $(s)$", fontsize=16) 
ax.set_ylabel(r"$g_2$", fontsize=16)



    


















    




















    
















    
Out[271]:








<matplotlib.text.Text at 0x7f9ff5f64748>

Another check for algorithm dealing with bad-frames





In [181]:



    


good_start = 6000
good_end = 8001   


good_start = 6420
good_end = 12801
#good_start = 7600
#good_end = 7650    

imgs_ =imgs[good_start: good_end-1]
imgsr = Reverse_Coordinate(imgs_, mask)



    











In [182]:



    


num_buf = 8



    











In [170]:



    


g2_t1, lag_steps_t1 = autocor_one_time( num_buf,  box_maskr, imgsr, num_lev=None,
                                bad_images=None, threshold= None )



    


















    






The lev number is 10
Doing g2 caculation of 2000 frames---
##########Total time: 5.05 min

















In [183]:



    


g2_t2, lag_steps_t2 = autocor_one_time( num_buf,  box_maskr, imgsr, num_lev=None,
                bad_images=None, threshold= 5000 )



    


















    






The lev number is 11
Doing g2 caculation of 6380 frames---
#bad image: 225 here!
#bad image: 825 here!
########bad image: 6164 here!
Total time: 16.18 min

















In [186]:



    


good_start = 8000
good_end = 12001
#good_start = 7600
#good_end = 7650    

imgs_ =imgs[good_start: good_end-1]
imgsr = Reverse_Coordinate(imgs_, mask)



    











In [187]:



    


g2_t4, lag_steps_t4 = autocor_one_time( num_buf,  box_maskr, imgsr, num_lev=None,
    bad_images=None, threshold= None )



    


















    






The lev number is 11
Doing g2 caculation of 4000 frames---
##########Total time: 10.28 min

















In [201]:



    


good_start = 7400
good_end = 12501
#good_start = 7600
#good_end = 7650    

imgs_ =imgs[good_start: good_end-1]
imgsr = Reverse_Coordinate(imgs_, mask)



    











In [202]:



    


g2_t5, lag_steps_t5 = autocor_one_time( num_buf,  box_maskr, imgsr, num_lev=None,
    bad_images=None, threshold= None )



    


















    






The lev number is 11
Doing g2 caculation of 5100 frames---
##########Total time: 12.94 min

















In [ ]:



    






    











In [188]:



    


g2_t2.shape



    


















    
Out[188]:








(45, 30)

















In [189]:



    


g2_t4.shape



    


















    
Out[189]:








(43, 30)

















In [197]:



    


plt.close('all')



    











In [206]:



    


fig,ax = plt.subplots()
ax.semilogx( lag_steps_t2*timeperframe, g2_t2[:,1], 'ro', label='inv-bad_images')
ax.semilogx( lag_steps_t5*timeperframe, g2_t5[:,1], '-b', label='all_good_images')
ax.set_ylim([ 1.1,1.2 ])
ax.legend()



    


















    




















    
















    
Out[206]:








<matplotlib.legend.Legend at 0x7f9fcea06828>

















In [ ]:



    






    











In [199]:



    


fig,ax = plt.subplots()
ax.semilogx( lag_steps_t2*timeperframe, g2_t2[:,1], 'ro', label='inv-bad_images')
ax.semilogx( lag_steps_t4*timeperframe, g2_t4[:,1], '-b', label='all_good_images')
ax.set_ylim([ 1.1,1.2 ])
ax.legend()



    


















    




















    
















    
Out[199]:








<matplotlib.legend.Legend at 0x7f9fceb06a90>

















In [ ]:



    






    











In [171]:



    


g2, lag_steps = g2_t1, lag_steps_t1
lags = lag_steps*timeperframe



    











In [184]:



    


g2, lag_steps = g2_t2, lag_steps_t2
lags = lag_steps*timeperframe



    











In [175]:



    


g2.shape



    


















    
Out[175]:








(39, 30)

















In [ ]:



    






    











In [172]:



    


tg2 = np.hstack( [ lags.reshape( len(lags),1), g2] )
#np.save( path + 'g2_%s-%s--%s'%(uid,good_start, good_end), tg2)
np.savetxt( path + 'g2_%s-%s--%s.txt'%(uid,good_start, good_end), tg2)

Plot the one time correlation functions





In [185]:



    


vmin = 1.1
vmax=  1.15

for qz_ind in range(num_qz):
    fig = plt.figure(figsize=(10, 12))
    #fig = plt.figure()
    title_qz = ' Qz= %.5f  '%( qz_center[qz_ind]) + r'$\AA^{-1}$' 
    plt.title('uid= %s:--->'%uid + title_qz,fontsize=20, y =1.1) 
    #print (qz_ind,title_qz)
    if num_qz!=1:plt.axis('off')
    sx = int(round(np.sqrt(num_qr)) )
    if num_qr%sx == 0: 
        sy = int(num_qr/sx)
    else:
        sy=int(num_qr/sx+1) 
    for sn in range(num_qr):
        ax = fig.add_subplot(sx,sy,sn+1 )
        ax.set_ylabel("g2") 
        ax.set_xlabel(r"$\tau $ $(s)$", fontsize=16) 
               
        title_qr = " Qr= " + '%.5f  '%( qr_center[sn]) + r'$\AA^{-1}$'
        if num_qz==1:
            
            title = 'uid= %s:--->'%uid + title_qz + '__' +  title_qr
        else:
            title = title_qr
        ax.set_title( title  )
            
            
        y=g2[:, sn + qz_ind * num_qr]
        ax.semilogx(lags, y, '-o', markersize=6) 
        ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
        #ax.set_ylim( vmin, vmax)
        
    plt.show()
    fig.tight_layout()  


#np.save( path + 'g2_%s--%s'%(good_start, good_end), [lags,g2])
#lags,g2 = np.load( path + 'g2_%s--%s.npy'%(good_start, good_end) )



    


















    






/home/yuzhang/.conda/envs/user_analysis/lib/python3.4/site-packages/matplotlib/pyplot.py:512: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  max_open_warning, RuntimeWarning)










    




















    
















    




















    























In [173]:



    


vmin = 1.1
vmax=  1.15

for qz_ind in range(num_qz):
    fig = plt.figure(figsize=(10, 12))
    #fig = plt.figure()
    title_qz = ' Qz= %.5f  '%( qz_center[qz_ind]) + r'$\AA^{-1}$' 
    plt.title('uid= %s:--->'%uid + title_qz,fontsize=20, y =1.1) 
    #print (qz_ind,title_qz)
    if num_qz!=1:plt.axis('off')
    sx = int(round(np.sqrt(num_qr)) )
    if num_qr%sx == 0: 
        sy = int(num_qr/sx)
    else:
        sy=int(num_qr/sx+1) 
    for sn in range(num_qr):
        ax = fig.add_subplot(sx,sy,sn+1 )
        ax.set_ylabel("g2") 
        ax.set_xlabel(r"$\tau $ $(s)$", fontsize=16) 
               
        title_qr = " Qr= " + '%.5f  '%( qr_center[sn]) + r'$\AA^{-1}$'
        if num_qz==1:
            
            title = 'uid= %s:--->'%uid + title_qz + '__' +  title_qr
        else:
            title = title_qr
        ax.set_title( title  )
            
            
        y=g2[:, sn + qz_ind * num_qr]
        ax.semilogx(lags, y, '-o', markersize=6) 
        ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
        #ax.set_ylim( vmin, vmax)
        
    plt.show()
    fig.tight_layout()  


#np.save( path + 'g2_%s--%s'%(good_start, good_end), [lags,g2])
#lags,g2 = np.load( path + 'g2_%s--%s.npy'%(good_start, good_end) )



    


















    




















    
















    




















    























In [174]:



    


vmin = 1.1
vmax=  1.15

for qz_ind in range(num_qz):
    fig = plt.figure(figsize=(10, 8))
    #fig = plt.figure()
    title_qz = ' Qz= %.5f  '%( qz_center[qz_ind]) + r'$\AA^{-1}$' 
    plt.title('uid= %s:--->'%uid + title_qz,fontsize=20, y =1.1) 
    #print (qz_ind,title_qz)
    if num_qz!=1:plt.axis('off')
    sx = int(round(np.sqrt(num_qr)) )
    if num_qr%sx == 0: 
        sy = int(num_qr/sx)
    else:
        sy=int(num_qr/sx+1) 
    for sn in range(num_qr):
        ax = fig.add_subplot(sx,sy,sn+1 )
        ax.set_ylabel("g2") 
        ax.set_xlabel(r"$\tau $ $(s)$", fontsize=16) 
               
        title_qr = " Qr= " + '%.5f  '%( qr_center[sn]) + r'$\AA^{-1}$'
        if num_qz==1:
            
            title = 'uid= %s:--->'%uid + title_qz + '__' +  title_qr
        else:
            title = title_qr
        ax.set_title( title  )
            
            
        y=g2[:, sn + qz_ind * num_qr]
        ax.semilogx(lags, y, '-o', markersize=6) 
        ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
        #ax.set_ylim( vmin, vmax)
        
    plt.show()
    fig.tight_layout()  


#np.save( path + 'g2_%s--%s'%(good_start, good_end), [lags,g2])
#lags,g2 = np.load( path + 'g2_%s--%s.npy'%(good_start, good_end) )

Fit g2





In [ ]:



    


from lmfit import  Model
mod = Model(corr.auto_corr_scat_factor)



    











In [ ]:



    


rate = []

sx = int( round (np.sqrt(num_rings)) )
if num_rings%sx==0:
    sy = int(num_rings/sx)
else:
    sy = int(num_rings/sx+1)
    
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid, fontsize=20, y =1.08)  
for sn in range(num_rings):
    ax = fig.add_subplot(sx, sy, sn+1 )
    y=g2[1:, sn]
    result1 = mod.fit(y, lags=lags[1:], beta=.1,
                      relaxation_rate =.5, baseline=1.0)
    rate.append(result1.best_values['relaxation_rate'])
    
    ax.semilogx(lags[1:], y, 'ro')
    ax.semilogx(lags[1:], result1.best_fit, '-b')
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')  
    ax.set_ylim([min(y)*.95, max(y[1:]) *1.05])
    txts = r'$\gamma$' + r'$ = %.3f$'%(rate[sn]) +  r'$ s^{-1}$'
    ax.text(x =0.015, y=.55, s=txts, fontsize=14, transform=ax.transAxes)              
 
fig.tight_layout()

Plot the relaxation rates vs (q_ring_center)**2





In [ ]:



    


fig, ax=plt.subplots()
ax.plot(q_ring_center**2, rate, 'ro', ls='--')
ax.set_ylabel('Relaxation rate 'r'$\gamma$'"($s^{-1}$)")
ax.set_xlabel("$q^2$"r'($\AA^{-2}$)')
plt.show()

Fitted the Diffusion Coefficinet D0





In [ ]:



    


D0 = np.polyfit(q_ring_center**2, rate, 1)
gmfit = np.poly1d(D0)
print ('The fitted diffusion coefficient D0 is:  %.2E   A^2S-1'%D0[0])



    











In [ ]:



    


fig,ax = plt.subplots()
ax.plot(q_ring_center**2, rate, 'ro', ls='')
ax.plot(q_ring_center**2,  gmfit(q_ring_center**2),  ls='-')
ax.set_ylabel('Relaxation rate 'r'$\gamma$'"($s^{-1}$)")
ax.set_xlabel("$q^2$"r'($\AA^{-2}$)')
plt.show()

Two_time Correlation





In [ ]:



    


good_start= 0
good_end = 10000
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)



    











In [ ]:



    


g12b = auto_two_Array( imgsr, ring_mask, data_pixel = data_pixel )



    











In [ ]:



    


g12_num = 0  #0: the firs ring
data = g12b[:,:,g12_num]
fig, ax = plt.subplots()
im=ax.imshow( data, origin='lower' , cmap='viridis', 
             norm= LogNorm( vmin= 1, vmax= 1.6 ), 
        extent=[0, data.shape[0]*timeperframe, 0, data.shape[0]*timeperframe ] )
ax.set_title('0-%s frames--Qth= %s'%(Nimg,g12_num))
ax.set_xlabel( r'$t_1$ $(s)$', fontsize = 18)
ax.set_ylabel( r'$t_2$ $(s)$', fontsize = 18)
fig.colorbar(im)
plt.show()

find the bad frames and mask it in two-time correlation





In [ ]:



    


#np.where( g12b[:,:,0] == g12b[:,:,0].min() )



    











In [ ]:



    


g12b_mask = make_g12_mask(g12b,[7629])
g12bm = masked_g12(g12b,[7629])

get one time correlation





In [ ]:



    


g2b = get_one_time_from_two_time(g12bm, timeperframe= timeperframe)



    











In [ ]:



    


sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0: 
    sy = int(num_rings/sx)
else:
    sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.08)  
plt.axis('off')
for sn in range(num_rings):
    ax = fig.add_subplot(sx,sy,sn+1 )
    ax.set_ylabel("g2") 
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')     
    y=g2b[:, sn]
    ax.semilogx( np.arange(len(y))*timeperframe, y, '-o', markersize=6) 
    #ax.semilogx(lags, y, '-o', markersize=6) 
    #ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
    ax.set_ylim( [1.0, max(y[1:])*1.05 ] )
    
plt.show()
fig.tight_layout()

plot g2 by multi-tau and g2 from two-time





In [ ]:



    


sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0: 
    sy = int(num_rings/sx)
else:
    sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.06)  
plt.axis('off')
for sn in range(num_rings):
    ax = fig.add_subplot(sx,sy,sn+1 )
    ax.set_ylabel("g2") 
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')     
    y=g2b[:, sn]
    ax.semilogx( np.arange(len(y))*timeperframe, y, '--r', markersize=6) 
    
    y2=g2[:, sn]
    ax.semilogx(lags, y2, 'o', markersize=6) 
    
    #ax.semilogx(lags, y, '-o', markersize=6) 
    #ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
    ax.set_ylim( [1.0, max(y2[1:])*1.05 ] )
plt.show()
fig.tight_layout()

Another defination of two-time correlation (different normalization)





In [ ]:



    


g12b_norm, g12b_not_norm, norms = auto_two_Array_g1_norm( imgsr, ring_mask, data_pixel = data_pixel )



    











In [ ]:



    


g2b_norm = get_one_time_from_two_time(g12b_norm,
                             norms=None, nopr=nopr, timeperframe= timeperframe)



    











In [ ]:



    


g2b_not_norm = get_one_time_from_two_time(g12b_not_norm,
                             norms=norms, nopr=nopr, timeperframe= timeperframe)



    











In [ ]:



    


#num_rings =6


sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0: 
    sy = int(num_rings/sx)
else:
    sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.06)  
plt.axis('off')
for sn in range(num_rings):
    ax = fig.add_subplot(sx,sy,sn+1 )
    ax.set_ylabel("g2") 
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')    
    
    y=g2b_norm[:, sn]
    #ax.semilogx( np.arange(len(y))*timeperframe, y, '--r', markersize=3) 
    
    y3=g2b_not_norm[:, sn]
    ax.semilogx( np.arange(len(y))*timeperframe, y3, '-g', lw=4, markersize=3)     
    
    y2=g2[:, sn]
    #ax.semilogx(lags, y2, '-ko', markersize= 3)
    ax.semilogx(lags, y2, 'o', markersize=6) 
    
    #ax.semilogx(lags, y, '-o', markersize=6) 
    #ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
    ax.set_ylim( [1.0, max(y2[1:])*1.05 ] )
plt.show()
fig.tight_layout()

get one time correlation @different age





In [ ]:



    


g12_num =0

g2b = show_g12_age_cuts( g12bm, g12_num=g12_num, slice_num =3, slice_width= 500, 
                slice_start=4000, slice_end= 20000-4000,
    timeperframe= timeperframe,vmin= 1.0, vmax =  1.22) #vmax=g12b[:,:,g12_num].max() )

get taus histgram





In [ ]:



    


taus = show_g12_tau_cuts(g12b_norm, g12_num=0,  slice_num = 5, slice_width= 1.0, 
    slice_start=3, slice_end= 5000-1,draw_scale_tau =1.0, vmin=1.01,vmax=1.55 )



    











In [ ]:



    


his = his_tau(taus, hisbin = 10, plot=True)



    











In [ ]:

Content source: yugangzhang/chx_backups

Similar notebooks:

notebook.community | gallery | about