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 [2]:

    

hdr = db[{{ uid }}]


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-27d71802241f> in <module>()
----> 1 hdr = db[{{ uid }}]

NameError: name 'uid' is not defined

In [3]:

    

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 [4]:

    

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


    

    




/home/yuzhang/chx-pipelines/Develops

In [ ]:

    

cd /XF11ID/analysis/Analysis_Pipelines/Develop/


    

In [7]:

    

%run develop.py
%run two_time.py


    

In [8]:

    

%matplotlib notebook
#%matplotlib inline


    

In [9]:

    

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


    

In [10]:

    

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 = 'ffe9d518'  # 10 mTorr 1sec/frame
else:
    uid = '/XF11ID/data/2015/11/23/d01ab510-3cf3-4719-bee3_795_master.h5'


    

In [14]:

    

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

In [12]:

    

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[12]:





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

In [15]:

    

imgs.md


    

    
Out[15]:





{'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 [16]:

    

# 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[16]:





(1.00001, 1.0)

In [17]:

    

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


    

In [18]:

    

maskr = mask[::-1,:]


    

In [19]:

    

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


    

In [20]:

    

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 [21]:

    

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


    

In [22]:

    

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)


    

In [23]:

    

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


    

In [24]:

    

bad_frames = np.where( kymo_sum > 1e5)[0]
bad_frames


    

    
Out[24]:





array([ 1124,  1245,  4304,  6645,  7245, 12584, 13784])

In [25]:

    

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


    

In [26]:

    

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


    

In [27]:

    

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


    

In [28]:

    

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()


    

In [29]:

    

# 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[29]:





(1.00001, 1.0)

In [30]:

    

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


    

In [31]:

    

lamda


    

    
Out[31]:





1.3776

In [32]:

    

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 [33]:

    

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 [34]:

    

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 [35]:

    

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(),) 
             #extent= [ qr.min(),qr.max(),qr.min(),qr.max()] )
fig.colorbar(im)
ax.set_title( 'Q-Parallel')
#ax.grid(True, which='both', color='r', linestyle='--')
ax.grid(True,  color='r', linestyle='--')
#ax.grid(True,color='white')

plt.show()


    

In [36]:

    

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()


    

In [37]:

    

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 [40]:

    

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


    

In [41]:

    

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


    

In [42]:

    

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


    

    















    












    















    












    




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

In [43]:

    

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


    

In [44]:

    

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 [45]:

    

num_qz,num_qr,num_qzr


    

    
Out[45]:





(2, 15, 30)

In [46]:

    

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 [47]:

    

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()


    

In [50]:

    

boxes = label_array_qzr 
box_maskr = boxes*maskr


    

In [51]:

    

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


    

In [52]:

    

nopr


    

    
Out[52]:





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 [56]:

    

# 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()


    

In [54]:

    

imgs_ =imgs
imgsr = Reverse_Coordinate(imgs_, maskr)


    

In [55]:

    

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


    

    




Total time: 4.40 min

In [57]:

    

max_inten_ring =2


    

In [58]:

    

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


    

In [59]:

    

data_pixel_qi.shape


    

    
Out[59]:





(15000, 15846)

In [60]:

    

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 [61]:

    

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()


    

In [62]:

    

mean_inten = get_mean_intensity( data_pixel, qind)


    

In [63]:

    

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 box")
for i in range(  num_qzr  ):
    ax.plot(times, mean_inten[i+1], '--o', label="Box "+str(i+1))
    ax.set_xlabel("Frame")
    ax.set_ylabel("Mean Intensity")
    ax.set_ylim(0, 2.5)
    
ax.legend( loc='best', fontsize = 6) 
plt.show()


    

In [64]:

    

mask_data =   np.arange( len(data_pixel))  
mask_data_ = np.delete( mask_data, bad_frames)


    

In [ ]:

    

#mask_data_


    

In [ ]:

    

#mean_inten[1][:1000].mean()


    

In [65]:

    

K_mean = np.array( [mean_inten[i][mask_data_].mean() for i in  list(mean_inten.keys() )] )


    

In [ ]:

    

K_mean


    

In [66]:

    

max_cts = data_pixel[mask_data_].max()
max_cts


    

    
Out[66]:





1180.0

In [67]:

    

data_pixel[1000].max()


    

    
Out[67]:





23.0

In [68]:

    

max_cts =60


    

In [69]:

    

data_pixel.shape


    

    
Out[69]:





(15000, 483222)

In [70]:

    

#max_cts = data_pixel.max()
#max_cts


    

In [71]:

    

%run speckle.py


    

In [72]:

    

good_start = 0
good_end =  3000# if using Nimg, then calculate all the images

good_end =    Nimg  #hen calculate all the images
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)


    

In [450]:

    

good_start = 4000
good_end =  8000# if using Nimg, then calculate all the images

#good_end =    Nimg  #hen calculate all the images
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)


    

In [451]:

    

len(imgs_)


    

    
Out[451]:





4000

In [452]:

    

time_steps = utils.geometric_series(2,   len(imgs_)   ) [:-3]
time_steps


    

    
Out[452]:





[1, 2, 4, 8, 16, 32, 64, 128, 256]

In [453]:

    

num_times = len(time_steps)


    

In [454]:

    

spe_cts_all, std_dev = xsvs( (imgsr,), np.int_(box_maskr), timebin_num=2,
        number_of_img= len(imgs_), max_cts=int(max_cts+2), time_bin = time_steps,
                            bad_images=None, threshold = 5000 )


    

    




#bad image: 304 here!
######bad image: 2645 here!
##bad image: 3245 here!
#Total time: 6.38 min

In [455]:

    

spe_cts_all.shape


    

    
Out[455]:





(9, 30)

In [ ]:

    

#np.save( path + 'uid_%s_spe_cts_all'%uid, spe_cts_all)


    

In [456]:

    

#np.save( path + 'uid_%s_spe_cts_%s--%s'%(uid,good_start,good_end), spe_cts_all)


    

In [76]:

    

#spe_cts_all = np.load(path + 'uid_%s_spe_cts_all.npy'%uid )


    

In [ ]:

    

#Knorm_bin_edges, Knorm_bin_centers = speckle.normalize_bin_edges( len(time_steps), num_rings, K_mean, int(max_cts+2))


    

In [457]:

    

bin_edges, bin_centers, Knorm_bin_edges, Knorm_bin_centers = get_bin_edges(
      len(time_steps), num_qzr, K_mean, int(max_cts+2)  )


    

In [458]:

    

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=16, y =1.05) 
    #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):
        axes = fig.add_subplot(sx,sy,sn+1 )        
        for j in range(len(time_steps)): 
            #if sn == 0:
            art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], '-o',
                     label=str(time_steps[j])+" ms")
            #else:
                #art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], '-o',)             
            
        axes.set_xlim(0, 3.5)         
        axes.legend(loc='best', fontsize = 6)
        
        axes.set_xlabel("K/<K>")
        axes.set_ylabel("P(K)")  
 
               
        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
        axes.set_title( title  )

    fig.tight_layout()


    

In [459]:

    

qz_ind =0
sn =3

fig = plt.figure(figsize=(8,6))
title_qz = ' Qz= %.5f  '%( qz_center[qz_ind]) + r'$\AA^{-1}$' 
#plt.title('uid= %s:--->'%uid + title_qz,fontsize=16, y =1.05) 
axes = fig.add_subplot(111 ) 

for j in range(len(time_steps)): 
    art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], '-o',
                     label=str(time_steps[j])+" ms")
            #else:
                #art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], '-o',)             
            
    axes.set_xlim(0, 5.0)         
    axes.legend(loc='best', fontsize = 12)
        
    axes.set_xlabel("K/<K>")
    axes.set_ylabel("P(K)")    

    title = 'uid= %s:--->'%uid + title_qz + '__' +  title_qr
    axes.set_title( title,fontsize=16, y =1.05) 

fig.tight_layout()


    

In [460]:

    

%run speckle.py


    

In [461]:

    

from lmfit import  Model
from scipy.interpolate import UnivariateSpline

#g_mod = Model(gamma_dist, indepdent_vars=['K'])
g_mod = Model( gamma_dist )
n_mod = Model(nbinom_dist)
p_mod = Model(poisson_dist)
dc_mod = Model(diff_mot_con_factor)


    

In [462]:

    

#gamma_dist??


    

In [463]:

    

#get_roi(data=spe_cts_all[0,3], threshold=1e-8)


    

In [466]:

    

qz_ind =0
sn = 11

fig = plt.figure(figsize=(8,6))
title_qz = ' Qz= %.5f  '%( qz_center[qz_ind]) + r'$\AA^{-1}$' 
#plt.title('uid= %s:--->'%uid + title_qz,fontsize=16, y =1.05) 
axes = fig.add_subplot(111 ) 

for j in range(len(time_steps)): 
    
    threshold=1e-4
    M= 10
    K= -1
    trys = 1
    for trys in range(7):
         
        if M<=0 or  K <= 0:
            threshold *= 10
            trys +=1
            
            n_mod.set_param_hint('M', value = 100, min=1.0, max= 10000 ) 
            n_mod.set_param_hint('K', value = K_mean[sn] * 2**j, 
                                 min=  K_mean[sn] * 2**j/10,max=  K_mean[sn] * 2**j*10) 
            
            #print (K_mean[sn]*2**j)
            
            
            pars = n_mod.make_params()                
            
            roi = get_roi(data=spe_cts_all[j, sn], threshold=threshold)      
            result_n = n_mod.fit(spe_cts_all[j, sn][roi],
                bin_values=bin_edges[j, sn][:-1][roi], params =pars) 
            M =  result_n.best_values['M']
            K =  result_n.best_values['K']
             
        else:
            break
            
           
    M_val=M
    K_val=K
    #print (M_val,K_val)
    
    #  Using the best K and M values interpolate and get more values for fitting curve
    #if j<=3:population_N = 10000
    #else:population_N=len(Knorm_bin_edges[j, sn][:-1])
    population_N  = 500*2**j   
    
    xlim=12.5
    fitx_ = np.linspace(0, max(Knorm_bin_edges[j, sn][:-1]), population_N     )   
    fitx = np.linspace(0, max(bin_edges[j, sn][:-1]), population_N )      
    
    #fitx_ = fitx_
    #fitx = np.linspace(0, xlim*K_mean[sn], population_N ) 
    fity = nbinom_dist( fitx, K_val, M_val ) # M and K are fitted best values
    
    
    if j == 0:
        art, = axes.plot( fitx_[:500],fity[:500], '-',  label="nbinom")
    else:
        art, = axes.plot( fitx_[:500],fity[:500], '-')                
                
    art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], 'o',
                     label=str(time_steps[j])+" ms")
            #else:
                #art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], '-o',)             
            
    axes.set_xlim(0, xlim)         
    axes.legend(loc='best', fontsize = 12)
   
    if j==0:axes.annotate(r'K='+'%.3f\n'%( K_val) +r'M='+'%.3f'%(M_val),
                      xy=(.25, 0.85),
                      xycoords='axes fraction', fontsize=14,
                      horizontalalignment='right', verticalalignment='bottom')
    axes.set_xlabel("K/<K>")
    axes.set_ylabel("P(K)")    

    title = 'uid= %s:--->'%uid + title_qz + '__' +  title_qr
    axes.set_title( title,fontsize=16, y =1.05) 

fig.tight_layout()


    

In [467]:

    

bn_M_val = {}
bn_K_val = {}

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=16, y =1.05) 
    #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) 
    #M_val[qz_ind]={}
    #K_val[qz_ind]={}
    for sn in range(num_qr):
        axes = fig.add_subplot(sx,sy,sn+1 ) 
        bn_M_val[qz_ind* num_qr + sn]=[]
        bn_K_val[qz_ind* num_qr + sn]=[]        
        for j in range(len(time_steps)):     
            threshold=1e-4
            M= 100
            K= -1
            trys = 1
            for trys in range(7):         
                if M<=0 or  K <= 0:
                    threshold *= 10
                    trys +=1
                    
                    K_ = K_mean[qz_ind* num_qr + sn]  * 2**j    
                    n_mod.set_param_hint('M', value = 100, min=1.0, max= 10000 ) 
                    n_mod.set_param_hint('K', value = K_, 
                                 min=  K_/10, max=  K_*10) 

                    pars = n_mod.make_params()     
                    roi = get_roi(data=spe_cts_all[j, sn], threshold=threshold)
                    
                    result_n = n_mod.fit(spe_cts_all[j, sn][roi],
                        bin_values=bin_edges[j, sn][:-1][roi], params =pars) 
                    M =  result_n.best_values['M']
                    K =  result_n.best_values['K']
             
                else:
                    break

            #print (M,K)
            bn_M_val[qz_ind* num_qr + sn].append(M)
            bn_K_val[qz_ind* num_qr + sn].append(K)
         
            #  Using the best K and M values interpolate and get more values for fitting curve
            population_N  = 500*2**j  
            #population_N  = 10000
            fitx_ = np.linspace(0, max(Knorm_bin_edges[j, sn][:-1]), population_N     )   
            fitx = np.linspace(0, max(bin_edges[j, sn][:-1]), population_N )  
            fity = nbinom_dist( fitx, bn_K_val[qz_ind* num_qr + sn][j],
                               bn_M_val[qz_ind* num_qr + sn][j] ) # M and K are fitted best values
 
            if j == 0:
                art, = axes.plot( fitx_[:500],fity[:500], '-',  label="nbinom")
            else:
                art, = axes.plot( fitx_[:500],fity[:500], '-') 
            #if j == 0:
               # art, = axes.plot( fitx_,fity, '-',  label="nbinom")
            #else:
                #art, = axes.plot( fitx_,fity, '-')        
                
            if j==0:axes.annotate(r'K='+'%.3f\n'%( K) +r'M='+'%.3f'%(M),
                      xy=(.65, 0.55),
                      xycoords='axes fraction', fontsize=14,
                      horizontalalignment='right', verticalalignment='bottom')
                
                
            if sn == 0:
                art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], 'o',
                     label=str(time_steps[j])+" ms")
            else:
                art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], 'o',)             
            
        axes.set_xlim(0, 3.5)         
        axes.legend(loc='best', fontsize = 6)
        
        axes.set_xlabel("K/<K>")
        axes.set_ylabel("P(K)")  
 
               
        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
        axes.set_title( title  )

    fig.tight_layout()


    

In [468]:

    

K_mean


    

    
Out[468]:





array([ 0.38237303,  0.17652289,  0.07161423,  0.04151951,  0.02951739,
        0.02423144,  0.02106524,  0.01891659,  0.01369414,  0.00836156,
        0.00419829,  0.00148836,  0.00123779,  0.00108712,  0.00100988,
        0.54096959,  0.18492724,  0.07727767,  0.04464761,  0.03106574,
        0.02452395,  0.02038103,  0.01715325,  0.0119164 ,  0.00718459,
        0.00355898,  0.00271693,  0.00226033,  0.00190928,  0.00168969])

In [469]:

    

K_mean_bn = []
for i in range(num_qz):
    for j in range(num_qr):
        K_mean_bn.append(    bn_K_val[i*num_qr+j][0]  )
K_mean_bn = np.array(K_mean_bn)


    

In [470]:

    

K_mean_bn


    

    
Out[470]:





array([ 0.37649687,  0.18294317,  0.07418293,  0.04289329,  0.03046884,
        0.02493529,  0.02163353,  0.01940351,  0.01403734,  0.00859343,
        0.00433231,  0.00155023,  0.00129475,  0.00113297,  0.00106373,
        0.37650529,  0.17996803,  0.07418293,  0.04289328,  0.03046884,
        0.02493529,  0.02163343,  0.01940351,  0.01403732,  0.00858802,
        0.00433231,  0.00155021,  0.00129475,  0.00113292,  0.00106372])

In [ ]:

    




    

In [472]:

    

qz_ind =0
sn =  3

fig = plt.figure(figsize=(8,6))
title_qz = ' Qz= %.5f  '%( qz_center[qz_ind]) + r'$\AA^{-1}$' 
#plt.title('uid= %s:--->'%uid + title_qz,fontsize=16, y =1.05) 
axes = fig.add_subplot(111 ) 

#for j in range(len(time_steps)): 
for j in range(4,len(time_steps)):   
    threshold=1e-6
    M= -10    
    trys = 1
    for trys in range(7):
         
        if M<=0:
            threshold *= 10
            trys +=1  
            M_ = bn_M_val[sn][j]
            #g_mod.set_param_hint('M', value = M_, min= M_*.01, max= M_*100 ) 
            g_mod.set_param_hint('M', value = 2.0 ) 
            #g_mod.set_param_hint('K', value = K_mean[sn] * 2**j, 
                                 #min=K_mean[sn] * 2**j*.99, max= K_mean[sn] * 2**j*1.01 ) 
            K_ = K_mean[sn]*2**j
            #print (M_,K_)
            g_mod.set_param_hint('K', value = K_, min=K_*.01, max= K_*100 )
            
            pars = g_mod.make_params()     
            roi = get_roi(data=spe_cts_all[j, sn], threshold=threshold)
                    
            result_g = g_mod.fit(spe_cts_all[j, sn][roi],
                    bin_values=bin_edges[j, sn][:-1][roi], params = pars)
                                 
                                 #parms = pars)
                            #K=K_mean[sn] * 2**j,
                                        # M =pars['M']) 
            M =  result_g.best_values['M'] 
            K =  result_g.best_values['K'] 
        else:
            break
            
           
    M_val=M
    K_val=K  #_mean[sn] * 2**j
    #print (M,K)
    #print (result_g.best_values['K'] )
    #print  (K_mean[sn] * 2**j)
            
    print (M_val,K_val)
    
    #  Using the best K and M values interpolate and get more values for fitting curve
    #if j<=3:population_N = 10000
    #else:population_N=len(Knorm_bin_edges[j, sn][:-1])
    population_N  = 5000*2**j   
    
    xlim=3.5
    fitx_ = np.linspace(0, max(Knorm_bin_edges[j, sn][:-1]), population_N     )   
    fitx = np.linspace(0, max(bin_edges[j, sn][:-1]), population_N )      
    
    #fitx_ = fitx_
    #fitx = np.linspace(0, xlim*K_mean[sn], population_N ) 
    fity = gamma_dist( fitx, K_val, M_val ) # M and K are fitted best values
    
    
    if j == 2:
        art, = axes.plot( fitx_[:5000],fity[:5000], '-k',  label="gamma")
    else:
        art, = axes.plot( fitx_[:5000],fity[:5000], '-r')                
                
    art, = axes.plot(Knorm_bin_edges[j, sn][:-1][roi], spe_cts_all[j, sn][roi], '--o',
                     label=str(time_steps[j])+" ms")
            #else:
                #art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], '-o',)             
            
    axes.set_xlim(0, 3)    
    axes.set_ylim(0, 0.3) 
    axes.legend(loc='best', fontsize = 12)
   
    if j==0:axes.annotate(r'K='+'%.3f\n'%( K_val) +r'M='+'%.3f'%(M_val),
                      xy=(.25, 0.85),
                      xycoords='axes fraction', fontsize=14,
                      horizontalalignment='right', verticalalignment='bottom')
    axes.set_xlabel("K/<K>")
    axes.set_ylabel("P(K)")    

    title = 'uid= %s:--->'%uid + title_qz + '__' +  title_qr
    axes.set_title( title,fontsize=16, y =1.05) 

fig.tight_layout()


    

    















    












    




0.999994242579 0.693292176881
1.3408048511 1.75448226949
2.72958892545 3.12665605525
4.65821094976 5.78568144535
8.53651588779 11.2642597124

In [365]:

    

%run speckle.py


    

In [473]:

    

gm_M_val = {}
gm_K_val={}

for qz_ind in range(num_qz):
    fig = plt.figure(figsize=(12,14))
    #fig = plt.figure()
    title_qz = ' Qz= %.5f  '%( qz_center[qz_ind]) + r'$\AA^{-1}$' 
    plt.title('uid= %s:--->'%uid + title_qz,fontsize=16, y =1.05) 
    #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) 
    #M_val[qz_ind]={}
    #K_val[qz_ind]={}
    
    for sn in range(num_qr):
        axes = fig.add_subplot(sx,sy,sn+1 ) 
        gm_M_val[qz_ind* num_qr + sn]=[]
        gm_K_val[qz_ind* num_qr + sn]=[]
        for j in range(len(time_steps)): 
        #for j in range(6,len(time_steps)):  
            threshold=1e-7
            M= -100            
            trys = 1
            for trys in range(7):         
                if M<=0:
                    threshold *= 10
                    trys +=1              
                    M_ = bn_M_val[sn][j]
                    K_ = K_mean[qz_ind* num_qr + sn] * 2**j
                    
                    g_mod.set_param_hint('M', value = 1.5)#, min= M_*.01, max= M_*100 )                     
                    #print (M_,K_)
                    g_mod.set_param_hint('K', value = K_, min=K_*.1, max= K_*10 )
 

                    pars = g_mod.make_params()     
                    roi = get_roi(data=spe_cts_all[j, sn], threshold=threshold)
                    
                    result_g = g_mod.fit(spe_cts_all[j, sn][roi],
                        bin_values=bin_edges[j, sn][:-1][roi],
                                        params = pars) 
                    M =  result_g.best_values['M']
                    K =  result_g.best_values['K']
                    #print (result_n.best_values['K']) 
             
                else:
                    break

            #print (M,K)
            gm_M_val[qz_ind* num_qr + sn].append(M)
            gm_K_val[qz_ind* num_qr + sn].append(K) 
         
            #  Using the best K and M values interpolate and get more values for fitting curve
            population_N  = 500*2**j  
            #population_N  = 10000
            fitx_ = np.linspace(0, max(Knorm_bin_edges[j, sn][:-1]), population_N     )   
            fitx = np.linspace(0, max(bin_edges[j, sn][:-1]), population_N )  
            fity = gamma_dist( fitx, gm_K_val[qz_ind* num_qr + sn][j], 
                              gm_M_val[qz_ind* num_qr + sn][j] ) # M and K are fitted best values
 
            if j == 0:
                art, = axes.plot( fitx_[:500],fity[:500], '-',  label="gamma")
            elif j>0:
                art, = axes.plot( fitx_[:500],fity[:500], '-') 
            else:
                pass
            #if j == 0:
               # art, = axes.plot( fitx_,fity, '-',  label="nbinom")
            #else:
                #art, = axes.plot( fitx_,fity, '-') 
        
            if sn == 0:
                art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], 'o',
                     label=str(time_steps[j])+" ms")
            else:
                art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], 'o',) 
                
            if j==0:
                axes.annotate(r'K='+'%.3f\n'%( K) +r'M='+'%.3f'%(M),
                      xy=(.65, 0.55),
                      xycoords='axes fraction', fontsize=14,
                      horizontalalignment='right', verticalalignment='bottom')
            
        axes.set_xlim(0, 3.5)      
        axes.set_ylim(0, 1)    
        axes.legend(loc='best', fontsize = 6)
        
        axes.set_xlabel("K/<K>")
        axes.set_ylabel("P(K)")  
 
               
        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
        axes.set_title( title  )

    fig.tight_layout()


    

In [300]:

    

#bin_edges[7,0]


    

In [474]:

    

qz_ind =0
sn =  3

fig = plt.figure(figsize=(8,6))
title_qz = ' Qz= %.5f  '%( qz_center[qz_ind]) + r'$\AA^{-1}$' 
#plt.title('uid= %s:--->'%uid + title_qz,fontsize=16, y =1.05) 
axes = fig.add_subplot(111 ) 

for j in range(len(time_steps)): 
#for j in range(6,len(time_steps)):   
    threshold=1e-8
    K= -10    
    trys = 1
    for trys in range(7):
         
        if K<=0:
            threshold *= 10
            trys +=1  
                        
            roi = get_roi(data=spe_cts_all[j, sn], threshold=threshold)
            K_ = K_mean[sn]*2**j
            p_mod.set_param_hint('K', value = K_, min=K_*.1, max= K_*10 )  
            pars = p_mod.make_params()    
            result_p = p_mod.fit(spe_cts_all[j, sn][roi],
                    bin_values=bin_edges[j, sn][:-1][roi],
                                 parms=pars)  
 
            K  =  result_p.best_values['K'] 
            #print (K,K_mean[sn]*2**j)
        else:
            break

    population_N  = 500*2**j   
    
    xlim=3.5
    fitx_ = np.linspace(0, max(Knorm_bin_edges[j, sn][:-1]), population_N     )   
    fitx = np.linspace(0, max(bin_edges[j, sn][:-1]), population_N )  
    
    #K=K_mean[sn]*2**j
    #print (K,K_mean[sn]*2**j)
    fity = poisson_dist( fitx, K  ) # M and K are fitted best values
    
    
    if j == 0:
        art, = axes.plot( fitx_[:500],fity[:500], '-r',  label="possion")
    else:
        art, = axes.plot( fitx_[:500],fity[:500], '-r')                
                
    #art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], 'o',
                     #label=str(time_steps[j])+" ms")
    art, = axes.plot(Knorm_bin_edges[j, sn][:-1][roi], spe_cts_all[j, sn][roi], 'o',
                     label=str(time_steps[j])+" ms")
            #else:
                #art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], '-o',)             
            
    axes.set_xlim(0, xlim)    
    #axes.set_ylim(0, 2) 
    axes.legend(loc='best', fontsize = 12)
   
    if j==0:axes.annotate(r'K='+'%.3f\n'%( K),
                      xy=(.25, 0.85),
                      xycoords='axes fraction', fontsize=14,
                      horizontalalignment='right', verticalalignment='bottom')
    axes.set_xlabel("K/<K>")
    axes.set_ylabel("P(K)")    

    title = 'uid= %s:--->'%uid + title_qz + '__' +  title_qr
    axes.set_title( title,fontsize=16, y =1.05) 

fig.tight_layout()


    

    















    












    




/home/yuzhang/.conda/envs/user_analysis/lib/python3.4/site-packages/lmfit-0.8.3-py3.4.egg/lmfit/model.py:455: UserWarning: The keyword argument parms does notmatch any arguments of the model function.It will be ignored.
  "It will be ignored.", UserWarning)

In [324]:

    

%run speckle.py


    

In [475]:

    

K_val = {}
for qz_ind in range(num_qz):
    fig = plt.figure(figsize=(12,14))
    #fig = plt.figure()
    title_qz = ' Qz= %.5f  '%( qz_center[qz_ind]) + r'$\AA^{-1}$' 
    plt.title('uid= %s:--->'%uid + title_qz,fontsize=16, y =1.05) 
    #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) 
    #M_val[qz_ind]={}
    #K_val[qz_ind]={}
    
    for sn in range(num_qr):
        axes = fig.add_subplot(sx,sy,sn+1 ) 
        K_val[qz_ind* num_qr + sn]=[]
         
        for j in range(len(time_steps)): 
        #for j in range(6,len(time_steps)):  
            threshold=1e-7
            K= -100            
            trys = 1
            for trys in range(7):         
                if K<=0:
                    threshold *= 10
                    trys +=1 
                    roi = get_roi(data=spe_cts_all[j, sn], threshold=threshold)
                    #K_ = K_mean[sn]*2**j
                    K_ = K_mean[qz_ind* num_qr + sn]  * 2**j   
                    p_mod.set_param_hint('K', value = K_, min=K_*.1, max= K_*10 )  
                    pars = p_mod.make_params()    
                    result_p = p_mod.fit(spe_cts_all[j, sn][roi],
                        bin_values=bin_edges[j, sn][:-1][roi],
                                 parms=pars)  
                    K  =  result_p.best_values['K']
                else:
                    break

            #print (M,K)
            K_val[qz_ind* num_qr + sn].append(K)
             
         
            #  Using the best K and M values interpolate and get more values for fitting curve
            population_N  = 500*2**j  
            #population_N  = 10000
            fitx_ = np.linspace(0, max(Knorm_bin_edges[j, sn][:-1]), population_N     )   
            fitx = np.linspace(0, max(bin_edges[j, sn][:-1]), population_N )  
            fity = poisson_dist( fitx,  K_val[qz_ind* num_qr + sn][j] ) # M and K are fitted best values

            if j == 0:
                art, = axes.plot( fitx_,fity, '-',  label="poisson")
            else:
                art, = axes.plot( fitx_,fity, '-') 
        
            if sn == 0:
                art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], 'o',
                     label=str(time_steps[j])+" ms")
            else:
                art, = axes.plot(Knorm_bin_edges[j, sn][:-1], spe_cts_all[j, sn], 'o',) 
                
            if j==0:axes.annotate(r'K='+'%.3f\n'%( K),
                      xy=(.65, 0.55),
                      xycoords='axes fraction', fontsize=14,
                      horizontalalignment='right', verticalalignment='bottom')
            
        axes.set_xlim(0, 3.5)         
        axes.legend(loc='best', fontsize = 6)
        
        axes.set_xlabel("K/<K>")
        axes.set_ylabel("P(K)")  

               
        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
        axes.set_title( title  )

    fig.tight_layout()


    

    















    












    




/home/yuzhang/.conda/envs/user_analysis/lib/python3.4/site-packages/lmfit-0.8.3-py3.4.egg/lmfit/model.py:455: UserWarning: The keyword argument parms does notmatch any arguments of the model function.It will be ignored.
  "It will be ignored.", UserWarning)






    















    

In [ ]:

    




    

In [476]:

    

contrast_factor = np.zeros((num_qzr, num_times))
for i in range(num_qz):
    for j in range(num_qr):
        contrast_factor[i*num_qr + j] =  1/np.array( bn_M_val[i*num_qr + j] )


    

In [477]:

    

contrast_factor[3  ]


    

    
Out[477]:





array([ 0.00516956,  0.01385574,  0.01665922,  0.02033683,  0.02526933,
        0.02312549,  0.0228302 ,  0.02373482,  0.02380551])

In [478]:

    

times = np.array( time_steps ) * exposuretime