"This notebook corresponds to version {{ version }} of the pipeline tool: https://github.com/NSLS-II/pipelines"
This notebook begins with a raw time-series of images and ends with $g_2(t)$ for a range of $q$, fit to an exponential or stretched exponential, and a two-time correlation functoin.
CHX Olog (https://logbook.nsls2.bnl.gov/11-ID/)
Import packages for I/O, visualization, and analysis.
In [2]:
from pyCHX.chx_packages import *
%matplotlib notebook
plt.rcParams.update({'figure.max_open_warning': 0})
plt.rcParams.update({ 'image.origin': 'lower' })
plt.rcParams.update({ 'image.interpolation': 'none' })
import pickle as cpk
from pyCHX.chx_xpcs_xsvs_jupyter_V1 import *
import itertools
from pyCHX.SAXS import poly_sphere_form_factor_intensity
In [3]:
#%matplotlib notebook
%matplotlib inline
In [4]:
scat_geometry = 'saxs' #suport 'saxs', 'gi_saxs', 'ang_saxs' (for anisotropics saxs or flow-xpcs)
#scat_geometry = 'ang_saxs'
#scat_geometry = 'gi_waxs'
#scat_geometry = 'gi_saxs'
roi_auto = True #False #True #if True, will automatically create a roi based on the roi_type ( iso/aniso etc), currently only works for SAXS
fit_peak = False #True
analysis_type_auto = True #False #True #False #True #if True, will take "analysis type" option from data acquisition func series
qphi_analysis = False #False #if True, will do q-phi (anisotropic analysis for transmission saxs)
isotropic_Q_mask = 'normal' #'wide' # 'normal' # 'wide' ## select wich Q-mask to use for rings: 'normal' or 'wide'
phi_Q_mask = 'phi_4x_20deg' ## select wich Q-mask to use for phi analysis
q_mask_name = ''
force_compress = False #True #force to compress data
bin_frame = False #generally make bin_frame as False
para_compress = True #parallel compress
run_fit_form = False #run fit form factor
run_waterfall = False #True #False #True #run waterfall analysis
run_profile_plot = False #run prolfile plot for gi-saxs
run_t_ROI_Inten = True #run ROI intensity as a function of time
run_get_mass_center = False # Analysis for mass center of reflective beam center
run_invariant_analysis = False
run_one_time = True #run one-time
cal_g2_error = False #True #calculate g2 signal to noise
#run_fit_g2 = True #run fit one-time, the default function is "stretched exponential"
fit_g2_func = 'stretched'
run_two_time = False #True #run two-time
run_four_time = False #True #True #False #run four-time
run_xsvs= False # True #False #run visibility analysis
att_pdf_report = True #attach the pdf report to CHX olog
qth_interest = 7 #the intested single qth
use_sqnorm = True #if True, use sq to normalize intensity
use_SG = True #False #True # False #if True, use the Sawitzky-Golay filter of the avg_img for normalization
use_SG_bin_frames = False #if True, use the Sawitzky-Golay filter of the (binned) frame for normalization
use_imgsum_norm= True #if True use imgsum to normalize intensity for one-time calculatoin
pdf_version='_%s'%get_today_date() #for pdf report name
run_dose = False # True # False #True # True #False #run dose_depend analysis
if scat_geometry == 'gi_saxs':run_xsvs= False;use_sqnorm=False
if scat_geometry == 'gi_waxs':use_sqnorm = False
if scat_geometry != 'saxs':qphi_analysis = False;scat_geometry_ = scat_geometry;roi_auto = False
else:scat_geometry_ = ['','ang_'][qphi_analysis]+ scat_geometry
if scat_geometry != 'gi_saxs':run_profile_plot = False
In [5]:
scat_geometry, scat_geometry_, qphi_analysis
Out[5]:
In [5]:
taus=None;g2=None;tausb=None;g2b=None;g12b=None;taus4=None;g4=None;times_xsv=None;contrast_factorL=None; lag_steps = None
In [6]:
CYCLE= '2019_3' #change clycle here
path = '/XF11ID/analysis/%s/masks/'%CYCLE
In [ ]:
In [9]:
username = getpass.getuser()
uid = '31930c82' #(scan num: 23309) (Measurement: test series for battery membrane )
print(uid)
username = 'koga'
roi_auto = True # True #False #True #if True, will automatically create a roi based on the roi_type ( iso/aniso etc), currently only works for SAXS
run_two_time = True
run_dose = False
In [10]:
#db[-1]
In [12]:
#get_last_uids( -1)
In [13]:
data_dir0 = create_user_folder(CYCLE, username)
print( data_dir0 )
In [14]:
uid = uid[:8]
print('The current uid for analysis is: %s...'%uid)
In [15]:
sud = get_sid_filenames(db[uid])
for pa in sud[2]:
if 'master.h5' in pa:
data_fullpath = pa
print ('scan_id, full-uid, data path are: %s--%s--%s'%(sud[0], sud[1], data_fullpath ))
#start_time, stop_time = '2017-2-24 12:23:00', '2017-2-24 13:42:00'
#sids, uids, fuids = find_uids(start_time, stop_time)
In [16]:
data_dir = os.path.join(data_dir0, '%s/'%(sud[1]))
os.makedirs(data_dir, exist_ok=True)
print('Results from this analysis will be stashed in the directory %s' % data_dir)
uidstr = 'uid=%s'%uid
In [17]:
md = get_meta_data( uid )
md_blue = md.copy()
#md_blue
In [18]:
#md_blue['detectors'][0]
#if md_blue['OAV_mode'] != 'none':
# cx , cy = md_blue[md_blue['detectors'][0]+'_beam_center_x'], md_blue[md_blue['detectors'][0]+'_beam_center_x']
#else:
# cx , cy = md_blue['beam_center_x'], md_blue['beam_center_y']
#print(cx,cy)
In [19]:
detectors = sorted(get_detectors(db[uid]))
print('The detectors are:%s'%detectors)
if len(detectors) >1:
md['detector'] = detectors[1]
print( md['detector'])
In [20]:
if md['detector'] =='eiger4m_single_image' or md['detector'] == 'image':
reverse= True
rot90= False
elif md['detector'] =='eiger500K_single_image':
reverse= True
rot90=True
elif md['detector'] =='eiger1m_single_image':
reverse= True
rot90=False
print('Image reverse: %s\nImage rotate 90: %s'%(reverse, rot90))
In [21]:
try:
cx , cy = md_blue['beam_center_x'], md_blue['beam_center_y']
print(cx,cy)
except:
print('Will find cx,cy later.')
In [22]:
if analysis_type_auto:#if True, will take "analysis type" option from data acquisition func series
try:
if scat_geometry != 'gi_saxs':
qphi_analysis_ = md['analysis'] #if True, will do q-phi (anisotropic analysis for transmission saxs)
print(md['analysis'])
if qphi_analysis_ == 'iso':
qphi_analysis = False
elif qphi_analysis_ == '':
qphi_analysis = False
else:
qphi_analysis = True
#for other analysis type, in case of GiSAXS, low_angle/high_anlge for instance
else:
gisaxs_inc_type = md['analysis']
except:
gisaxs_inc_type = None
print('There is no analysis in metadata.')
print('Will %s qphis analysis.'%['NOT DO','DO'][qphi_analysis])
if scat_geometry != 'saxs':qphi_analysis = False;scat_geometry_ = scat_geometry
else:scat_geometry_ = ['','ang_'][qphi_analysis]+ scat_geometry
if scat_geometry != 'gi_saxs':run_profile_plot = False
print(scat_geometry_)
In [23]:
scat_geometry
Out[23]:
In [24]:
#%run /home/yuzhang/pyCHX_link/pyCHX/chx_generic_functions.py
In [25]:
imgs = load_data( uid, md['detector'], reverse= reverse, rot90=rot90 )
md.update( imgs.md );Nimg = len(imgs);
#md['beam_center_x'], md['beam_center_y'] = cx, cy
#if 'number of images' not in list(md.keys()):
md['number of images'] = Nimg
pixel_mask = 1- np.int_( np.array( imgs.md['pixel_mask'], dtype= bool) )
print( 'The data are: %s' %imgs )
#md['acquire period' ] = md['cam_acquire_period']
#md['exposure time'] = md['cam_acquire_time']
mdn = md.copy()
In [26]:
if md['detector'] =='eiger1m_single_image':
Chip_Mask=np.load( '/XF11ID/analysis/2017_1/masks/Eiger1M_Chip_Mask.npy')
elif md['detector'] =='eiger4m_single_image' or md['detector'] == 'image':
Chip_Mask= np.array(np.load( '/XF11ID/analysis/2017_1/masks/Eiger4M_chip_mask.npy'), dtype=bool)
BadPix = np.load('/XF11ID/analysis/2018_1/BadPix_4M.npy' )
Chip_Mask.ravel()[BadPix] = 0
elif md['detector'] =='eiger500K_single_image':
#print('here')
Chip_Mask= np.load( '/XF11ID/analysis/2017_1/masks/Eiger500K_Chip_Mask.npy') #to be defined the chip mask
Chip_Mask = np.rot90(Chip_Mask)
pixel_mask = np.rot90( 1- np.int_( np.array( imgs.md['pixel_mask'], dtype= bool)) )
else:
Chip_Mask = 1
print(Chip_Mask.shape, pixel_mask.shape)
In [27]:
use_local_disk = True
import shutil,glob
In [28]:
save_oavs = False
if len(detectors)==2:
if '_image' in md['detector']:
pref = md['detector'][:-5]
else:
pref=md['detector']
for k in [ 'beam_center_x', 'beam_center_y','cam_acquire_time','cam_acquire_period','cam_num_images',
'wavelength', 'det_distance', 'photon_energy']:
md[k] = md[ pref + '%s'%k]
if 'OAV_image' in detectors:
try:
save_oavs_tifs( uid, data_dir )
save_oavs = True
except:
pass
In [29]:
print_dict( md, ['suid', 'number of images', 'uid', 'scan_id', 'start_time', 'stop_time', 'sample', 'Measurement',
'acquire period', 'exposure time',
'det_distance', 'beam_center_x', 'beam_center_y', ] )
In [30]:
if scat_geometry =='gi_saxs':
inc_x0 = md['beam_center_x']
inc_y0 = imgs[0].shape[0] - md['beam_center_y']
#1147 ,1550
refl_x0 = md['beam_center_x']
refl_y0 = 1550 # imgs[0].shape[0] - 673
print( "inc_x0, inc_y0, ref_x0,ref_y0 are: %s %s %s %s."%(inc_x0, inc_y0, refl_x0, refl_y0) )
else:
if md['detector'] =='eiger4m_single_image' or md['detector'] == 'image' or md['detector']=='eiger1m_single_image':
inc_x0 = imgs[0].shape[0] - md['beam_center_y']
inc_y0= md['beam_center_x']
elif md['detector'] =='eiger500K_single_image':
inc_y0 = imgs[0].shape[1] - md['beam_center_y']
inc_x0 = imgs[0].shape[0] - md['beam_center_x']
print(inc_x0, inc_y0)
###for this particular uid, manually give x0/y0
#inc_x0 = 1041
#inc_y0 = 1085
In [31]:
dpix, lambda_, Ldet, exposuretime, timeperframe, center = check_lost_metadata(
md, Nimg, inc_x0 = inc_x0, inc_y0= inc_y0, pixelsize = 7.5*10*(-5) )
if scat_geometry =='gi_saxs':center=center[::-1]
if uid == '4740b6b5-20cb-4ad1-9f7a-b7a2b3a77659':
Ldet = 16.035*1000
setup_pargs=dict(uid=uidstr, dpix= dpix, Ldet=Ldet, lambda_= lambda_, exposuretime=exposuretime,
timeperframe=timeperframe, center=center, path= data_dir)
print_dict( setup_pargs )
In [32]:
setup_pargs
Out[32]:
In [33]:
date = 'Sep30'
if scat_geometry == 'gi_saxs':
mask_path = '/XF11ID/analysis/2019_3/masks/'
mask_name = '%s_2019_4M_SAXS.npy'%date
elif scat_geometry == 'saxs':
mask_path = '/XF11ID/analysis/2019_3/masks/'
if md['detector'] =='eiger4m_single_image' or md['detector'] == 'image':
mask_name = '%s_2019_4M_SAXS.npy'%date
elif md['detector'] =='eiger500K_single_image':
mask_name = '%s_2019_500K_SAXS.npy'%date
elif scat_geometry == 'gi_waxs':
mask_path = '/XF11ID/analysis/2019_3/masks/'
mask_name = '%s_2019_1M_WAXS.npy'
In [ ]:
In [34]:
mask = load_mask(mask_path, mask_name, plot_ = False, image_name = uidstr + '_mask', reverse= reverse, rot90=rot90 )
mask = mask * pixel_mask * Chip_Mask
show_img(mask,image_name = uidstr + '_mask', save=True, path=data_dir, aspect=1, center=center[::-1])
mask_load=mask.copy()
imgsa = apply_mask( imgs, mask )
In [35]:
#%run ~/pyCHX_link/pyCHX/chx_generic_functions.py
In [36]:
print(roi_auto, qphi_analysis, isotropic_Q_mask)
In [37]:
roi_date = 'Sep13'
if scat_geometry =='saxs':
## For auto load roi mask
if roi_auto:
general_path = '/nsls2/xf11id1/analysis/masks/'
roi_mask_norm = general_path + 'general_roi_mask_norm.npy'
roi_mask_wide = general_path + 'general_roi_mask_wide.npy'
roi_mask_phi_4x_20deg = general_path + 'general_roi_mask_phi_4x_20deg.npy'
roi_mask_12x_30deg = general_path + 'general_roi_mask_12x_30deg.npy'
roi_mask_12x_15deg_flow = general_path + 'general_roi_mask_12x_15deg_flow.npy'
if qphi_analysis == False:
if isotropic_Q_mask == 'normal':
fp = roi_mask_norm
elif isotropic_Q_mask == 'wide':
fp = roi_mask_wide
elif qphi_analysis:
if phi_Q_mask == 'phi_4x_20deg':
fp = roi_mask_phi_4x_20deg
elif phi_Q_mask == 'phi_12x_30deg':
fp = roi_mask_phi_12x_30deg
elif phi_Q_mask == 'phi_12x_15deg_flow':
fp = roi_mask_phi_12x_15deg_flow
roi_mask0 = np.load(fp)
old_cen=[4000,4000]
roi_mask, qval_dict, qwid_dict = get_roi_mask_qval_qwid_by_shift(
new_cen=center, new_mask= mask, old_cen=old_cen,
old_roi_mask=roi_mask0, limit_qnum= None,
setup_pargs= setup_pargs, geometry = scat_geometry_,
)
else:
if qphi_analysis == False:
if isotropic_Q_mask == 'normal':
#print('Here')
q_mask_name='rings'
if md['detector'] =='eiger4m_single_image' or md['detector'] == 'image': #for 4M
fp = roi_path + 'roi_mask_%s_4M_norm.pkl'%roi_date
elif md['detector'] =='eiger500K_single_image': #for 500K
fp = roi_path + 'roi_mask_%s_500K_norm.pkl'%roi_date
elif isotropic_Q_mask == 'wide':
q_mask_name='wide_rings'
if md['detector'] =='eiger4m_single_image' or md['detector'] == 'image': #for 4M
fp = roi_path + 'roi_mask_%s_4M_wide.pkl'%roi_date
elif md['detector'] =='eiger500K_single_image': #for 500K
fp = roi_path + 'roi_mask_%s_500K_wide.pkl'%roi_date
elif qphi_analysis:
if phi_Q_mask =='phi_4x_20deg':
q_mask_name='phi_4x_20deg'
if md['detector'] =='eiger4m_single_image' or md['detector'] == 'image': #for 4M
fp = roi_path + 'roi_mask_%s_4M_phi_4x_20deg.pkl'%roi_date
fp = roi_path + 'roi_mask_%s_4M_phi_24x_15deg.pkl'%roi_date
elif md['detector'] =='eiger500K_single_image': #for 500K
fp = roi_path + 'roi_mask_%s_500K_phi_4x_20deg.pkl'%roi_date
#fp = 'XXXXXXX.pkl'
roi_mask,qval_dict = cpk.load( open(fp, 'rb' ) ) #for load the saved roi data
print(fp)
## Gi_SAXS
elif scat_geometry =='gi_saxs':
# static mask
ss = '3e4f076a'
fp = '/XF11ID/analysis/2019_3/masks/uid=%s_roi_masks.pkl'%ss
roi_masks,qval_dicts = cpk.load( open(fp, 'rb' ) ) #for load the saved roi data
print('The static mask is: %s.'%fp)
fp = '/XF11ID/analysis/2019_3/masks/uid=%s_roi_mask.pkl'%(ss)
roi_mask,qval_dict = cpk.load( open(fp, 'rb' ) ) #for load the saved roi data
print('The dynamic mask is: %s.'%fp)
# q-map
fp = '/XF11ID/analysis/2019_3/masks/uid=%s_qmap.pkl'%(ss)
#print(fp)
qr_map, qz_map, ticks, Qrs, Qzs, Qr, Qz, inc_x0,refl_x0, refl_y0 = cpk.load( open(fp, 'rb' ) )
print('The qmap is: %s.'%fp)
## WAXS
elif scat_geometry =='gi_waxs':
#fp = '/XF11ID/analysis/2019_1/masks/uid=ca2ccb14_roi_mask_5PTO_130C_PTO.pkl'
fp = '/XF11ID/analysis/2019_2/masks/uid=278c0df1_roi_mask.pkl'
roi_mask,qval_dict = cpk.load( open(fp, 'rb' ) ) #for load the saved roi data
print(roi_mask.shape)
In [38]:
#qval_dict
In [39]:
#roi_mask = shift_mask(roi_mask, 10,30) #if shift mask to get new mask
show_img(roi_mask, aspect=1.0, image_name = fp)#, center=center[::-1])
In [40]:
img_choice_N = 1
img_samp_index = random.sample( range(len(imgs)), img_choice_N)
avg_img = get_avg_img( imgsa, img_samp_index, plot_ = False, uid =uidstr)
if avg_img.max() == 0:
print('There are no photons recorded for this uid: %s'%uid)
print('The data analysis should be terminated! Please try another uid.')
In [41]:
#show_img( imgsa[1000], vmin=.1, vmax= 1e1, logs=True, aspect=1,
# image_name= uidstr + '_img_avg', save=True, path=data_dir, cmap = cmap_albula )
print(center[::-1])
In [42]:
show_img( imgsa[ 5], vmin = 1, vmax = 10, logs=False, aspect=1, #save_format='tif',
image_name= uidstr + '_img_avg', save=True, path=data_dir, cmap=cmap_albula,center=center[::-1])
# select subregion, hard coded center beam location
#show_img( imgsa[180+40*3/0.05][110:110+840*2, 370:370+840*2], vmin = 0.01, vmax = 20, logs=False, aspect=1, #save_format='tif',
# image_name= uidstr + '_img_avg', save=True, path=data_dir, cmap=cmap_albula,center=[845,839])
In [43]:
compress=True
photon_occ = len( np.where(avg_img)[0] ) / ( imgsa[0].size)
#compress = photon_occ < .4 #if the photon ocupation < 0.5, do compress
print ("The non-zeros photon occupation is %s."%( photon_occ))
print("Will " + 'Always ' + ['NOT', 'DO'][compress] + " apply compress process.")
In [44]:
if md['detector'] =='eiger4m_single_image' or md['detector'] == 'image':
good_start = 5 #make the good_start at least 0
elif md['detector'] =='eiger500K_single_image':
good_start = 5#100 #5 #make the good_start at least 0
elif md['detector'] =='eiger1m_single_image' or md['detector'] == 'image':
good_start = 5
In [45]:
bin_frame = False # True #generally make bin_frame as False
if bin_frame:
bin_frame_number=4
acquisition_period = md['acquire period']
timeperframe = acquisition_period * bin_frame_number
else:
bin_frame_number =1
In [46]:
force_compress = False
#force_compress = True
#force_compress = True
In [47]:
import time
t0= time.time()
if not use_local_disk:
cmp_path = '/nsls2/xf11id1/analysis/Compressed_Data'
else:
cmp_path = '/tmp_data/compressed'
cmp_path = '/nsls2/xf11id1/analysis/Compressed_Data'
if bin_frame_number==1:
cmp_file = '/uid_%s.cmp'%md['uid']
else:
cmp_file = '/uid_%s_bined--%s.cmp'%(md['uid'],bin_frame_number)
filename = cmp_path + cmp_file
mask2, avg_img, imgsum, bad_frame_list = compress_eigerdata(imgs, mask, md, filename,
force_compress= force_compress, para_compress= para_compress, bad_pixel_threshold = 1e14,
reverse=reverse, rot90=rot90, nobytes=4,
bins=bin_frame_number, num_sub= 100, num_max_para_process= 500, with_pickle=True,
direct_load_data =use_local_disk, data_path = data_fullpath, )
min_inten = 10
good_start = max(good_start, np.where( np.array(imgsum) > min_inten )[0][0] )
print ('The good_start frame number is: %s '%good_start)
FD = Multifile(filename, good_start, len(imgs)//bin_frame_number )
#FD = MultifileBNLCustom(filename, good_start, len(imgs)//bin_frame_number )
###For test purpose to use the first 1000 frames
####################################################################
#FD = MultifileBNLCustom(filename, good_start, 500 )
#FD = Multifile(filename, good_start, 500)
############################################################
uid_ = uidstr + '_fra_%s_%s'%(FD.beg, FD.end)
print( uid_ )
plot1D( y = imgsum[ np.array( [i for i in np.arange(good_start, len(imgsum)) if i not in bad_frame_list])],
title =uidstr + '_imgsum', xlabel='Frame', ylabel='Total_Intensity', legend='imgsum' )
Nimg = int( Nimg/bin_frame_number )
run_time(t0)
mask2 = mask * pixel_mask * Chip_Mask
mask_copy = mask.copy()
mask_copy2 = mask.copy()
avg_img *= mask
In [48]:
#%run ~/pyCHX_link/pyCHX/chx_generic_functions.py
In [49]:
try:
if md['experiment']=='printing':
#p = md['printing'] #if have this printing key, will do error function fitting to find t_print0
find_tp0 = True
t_print0 = ps( y = imgsum[:400] ) * timeperframe
print( 'The start time of print: %s.' %(t_print0 ) )
else:
find_tp0 = False
print('md[experiment] is not "printing" -> not going to look for t_0')
t_print0 = None
except:
find_tp0 = False
print('md[experiment] is not "printing" -> not going to look for t_0')
t_print0 = None
In [50]:
show_img( avg_img*mask, vmin=1e-5, vmax= 1e3, logs=True, aspect=1, #save_format='tif',
image_name= uidstr + '_img_avg', save=True,
path=data_dir, center=center[::-1], cmap = cmap_albula ) #plt.cm.jet )#
In [51]:
good_end= None # 2000
if good_end is not None:
FD = Multifile(filename, good_start, min( len(imgs)//bin_frame_number, good_end) )
uid_ = uidstr + '_fra_%s_%s'%(FD.beg, FD.end)
print( uid_ )
In [52]:
re_define_good_start =False
if re_define_good_start:
good_start = 180
#good_end = 19700
good_end = len(imgs)
FD = Multifile(filename, good_start, good_end)
uid_ = uidstr + '_fra_%s_%s'%(FD.beg, FD.end)
print( FD.beg, FD.end)
In [53]:
bad_frame_list = get_bad_frame_list( imgsum, fit='both', plot=True,polyfit_order = 30,
scale= 3.5, good_start = good_start, good_end=good_end, uid= uidstr, path=data_dir)
print( 'The bad frame list length is: %s'%len(bad_frame_list) )
In [54]:
imgsum_y = imgsum[ np.array( [i for i in np.arange( len(imgsum)) if i not in bad_frame_list])]
imgsum_x = np.arange( len( imgsum_y))
save_lists( [imgsum_x, imgsum_y], label=['Frame', 'Total_Intensity'],
filename=uidstr + '_img_sum_t', path= data_dir )
In [55]:
plot1D( y = imgsum_y, title = uidstr + '_img_sum_t', xlabel='Frame', c='b',
ylabel='Total_Intensity', legend='imgsum', save=True, path=data_dir)
In [56]:
if md['detector'] =='eiger4m_single_image' or md['detector'] == 'image':
pass
elif md['detector'] =='eiger500K_single_image':
#if md['cam_acquire_period'] <= 0.00015: #will check this logic
if imgs[0].dtype == 'uint16':
print('Create dynamic mask for 500K due to 9K data acquistion!!!')
#bdp = find_bad_pixels_FD( bad_frame_list, FD, img_shape = avg_img.shape, threshold= 15 )
bdp = find_bad_pixels_FD( np.arange(Nimg), FD, img_shape = avg_img.shape, threshold= 15 )
mask = mask_copy2.copy()
mask *=bdp
mask_copy = mask.copy()
show_img( mask, image_name='New Mask_uid=%s'%uid, aspect=1 )
In [57]:
#show_img( mask )
In [ ]:
In [58]:
setup_pargs
Out[58]:
In [59]:
if scat_geometry =='saxs':
## Get circular average| * Do plot and save q~iq
mask = mask_copy.copy()
hmask = create_hot_pixel_mask( avg_img, threshold = 1e8, center=center, center_radius= 10)
qp_saxs, iq_saxs, q_saxs = get_circular_average( avg_img * Chip_Mask , mask * hmask, save=True,
pargs=setup_pargs )
plot_circular_average( qp_saxs, iq_saxs, q_saxs, pargs=setup_pargs, show_pixel=True,
xlim=[qp_saxs.min(), qp_saxs.max()*1.0], ylim = [iq_saxs.min(), iq_saxs.max()*2] )
mask =np.array( mask * hmask, dtype=bool)
In [60]:
#%run -i ~/pyCHX_link/pyCHX/SAXS.py
In [61]:
if scat_geometry =='saxs':
if run_fit_form:
qp, pq = q_saxs, iq_saxs
fitp, fitpe, q_, fitq, result = get_form_factor_fit_lmfit( qp, pq,
guess_values={'radius': 1180, 'sigma':0.05,'delta_rho':1e-9, 'background':5e-5,
'qpower_scale':5e-10, 'qpower': -4.0 },
guess_limit = { 'radius': [1000, 1500], 'sigma':[1e-3, .5], 'delta_rho': [1e-16, 1e-4],
'background': [1e-6, 1e-3], 'qpower_scale': [1e-16, 1], 'qpower': [ -5,-3]},
fit_range=[0.005, 0.03],
fit_variables={'radius': True, 'sigma':True, 'delta_rho':True, 'background': True,
'qpower_scale':True, 'qpower':False},
fit_func='G', qpower_bg= True )
r,sig=fitp[ 'radius'], fitp['sigma']
radius, sigma, delta_rho, background, qpower_scale, qpower = (fitp[ 'radius'], fitp['sigma'],
fitp['delta_rho'],fitp['background'],
fitp[ 'qpower_scale'], fitp['qpower'] )
bk= qpower_scale*qp**( qpower ) + background
ff= poly_sphere_form_factor_intensity( qp, radius=radius, sigma=sigma, delta_rho=delta_rho,
background=0, num_points= 20,
spread=5,fit_func= 'G' ) + bk
#report_fit( result)
In [62]:
if scat_geometry =='saxs':
if run_fit_form:
fig = plt.figure( )#figsize=[8,6])
ax = fig.add_subplot( 1,1,1 )
plot1D(x=qp,y=pq,ax=ax,m='o', ls='',c='k', legend='data', logy=True)
#plot1D(x=q_,y=fitq,ax=ax,m='', ls='-',c='r', legend='fitG_bg', logy=True)#, logx=True )
plot1D(x=qp,y=ff,ax=ax,m='', ls='-',c='b', legend='fitG_bgA', logy=True)#, logx=True )
#plot1D(x=q, y = bk, ls='--',c='g',ax=ax, lw=1, markersize=2,)
#qpower_scale, qpower
plot1D(x=qp, y= qp**( qpower ) * np.exp( qpower_scale ) /1e8, ax=ax, m='', ls='-', c='k', legend= r'$q^{%.1f}$'% qpower)
txts = r'radius' + r' = %.2f +/- %.2f '%( radius/10, fitpe['radius'] ) + r'$ nm$'
ax.text(x =0.02, y=.2, s=txts, fontsize=14, transform=ax.transAxes)
txts = r'sigma' + r' = %.4f +/- %.4f'%( sigma, fitpe['sigma'])
#txts = r'$\beta$' + r'$ = %.3f$'%(beta[i]) + r'$ s^{-1}$'
ax.text(x =0.02, y=.1, s=txts, fontsize=14, transform=ax.transAxes)
ax.set_xscale('log')
#ax.set_ylim( [1.00, 1.15 ] )
ax.set_xlabel(r'$q$ ('r'$\AA^{-1}$)', fontsize=18)
ax.set_ylabel(r'$I(q)$', fontsize=18)
#ax.set_title( 'uid=%s'%uid, fontsize=18)
plt.tight_layout()
In [63]:
#qval_dict
In [64]:
if scat_geometry =='saxs':
qr = np.array( [qval_dict[k][0] for k in sorted( qval_dict.keys())] )
if qphi_analysis == False:
try:
qr_cal, qr_wid = get_QrQw_From_RoiMask( roi_mask, setup_pargs )
print(len(qr))
if (qr_cal - qr).sum() >=1e-3:
print( 'The loaded ROI mask might not be applicable to this UID: %s.'%uid)
print('Please check the loaded roi mask file.')
except:
print('Something is wrong with the roi-mask. Please check the loaded roi mask file.')
show_ROI_on_image( avg_img*roi_mask, roi_mask, center, label_on = False, rwidth = 840, alpha=.9,
save=True, path=data_dir, uid=uidstr, vmin= 1e-3,
vmax= 1e-1, #np.max(avg_img),
aspect=1,
show_roi_edge=True,
show_ang_cor = True)
plot_qIq_with_ROI( q_saxs, iq_saxs, np.unique(qr), logs=True, uid=uidstr,
xlim=[q_saxs.min(), q_saxs.max()*1.02],#[0.0001,0.08],
ylim = [iq_saxs.min(), iq_saxs.max()*1.02], save=True, path=data_dir)
roi_mask = roi_mask * mask
In [ ]:
In [65]:
if scat_geometry =='saxs':
Nimg = FD.end - FD.beg
time_edge = create_time_slice( Nimg, slice_num= 3, slice_width= 2, edges = None )
time_edge = np.array( time_edge ) + good_start
#print( time_edge )
qpt, iqst, qt = get_t_iqc( FD, time_edge, mask*Chip_Mask, pargs=setup_pargs, nx=1500, show_progress= False )
plot_t_iqc( qt, iqst, time_edge, pargs=setup_pargs, xlim=[qt.min(), qt.max()],
ylim = [iqst.min(), iqst.max()], save=True )
In [66]:
#fit_peak = True
#from pyCHX.chx_generic_functions import fit_one_peak_curve
#%run -i ~/pyCHX_link/pyCHX/chx_generic_functions.py
In [67]:
if scat_geometry =='saxs':
if fit_peak:
q1 = 0.162
x = q_saxs
y= iq_saxs
aind = [ find_index( q_saxs, q1 ) ]
w = 50
xrange = [ aind[0] -w, aind[0] + w]
cen, cen_std, wid, wid_std, xf, out = fit_one_peak_curve( x, y, xrange )
fp = data_dir + uidstr + '_qIq_fit.png'
ax = plot_xy_with_fit( x, y, xf, out, cen, cen_std,wid, wid_std,
xlim=[0.01,0.3],xlabel= 'q ('r'$\AA^{-1}$)',
ylabel='I(q)', filename=fp)
In [68]:
if run_invariant_analysis:
if scat_geometry =='saxs':
invariant = get_iq_invariant( qt, iqst )
time_stamp = time_edge[:,0] * timeperframe
if scat_geometry =='saxs':
plot_q2_iq( qt, iqst, time_stamp,pargs=setup_pargs,ylim=[ -0.001, 0.01] ,
xlim=[0.007,0.2],legend_size= 6 )
if scat_geometry =='saxs':
plot_time_iq_invariant( time_stamp, invariant, pargs=setup_pargs, )
if False:
iq_int = np.zeros( len(iqst) )
fig, ax = plt.subplots()
q = qt
for i in range(iqst.shape[0]):
yi = iqst[i] * q**2
iq_int[i] = yi.sum()
time_labeli = 'time_%s s'%( round( time_edge[i][0] * timeperframe, 3) )
plot1D( x = q, y = yi, legend= time_labeli, xlabel='Q (A-1)', ylabel='I(q)*Q^2', title='I(q)*Q^2 ~ time',
m=markers[i], c = colors[i], ax=ax, ylim=[ -0.001, 0.01] , xlim=[0.007,0.2],
legend_size=4)
#print( iq_int )
In [69]:
if scat_geometry =='gi_saxs':
plot_qzr_map( qr_map, qz_map, inc_x0, ticks = ticks, data= avg_img, uid= uidstr, path = data_dir )
In [70]:
if scat_geometry =='gi_saxs':
#roi_masks, qval_dicts = get_gisaxs_roi( Qrs, Qzs, qr_map, qz_map, mask= mask )
show_qzr_roi( avg_img, roi_masks, inc_x0, ticks[:4], alpha=0.5, save=True, path=data_dir, uid=uidstr )
In [71]:
if scat_geometry =='gi_saxs':
Nimg = FD.end - FD.beg
time_edge = create_time_slice( N= Nimg, slice_num= 3, slice_width= 2, edges = None )
time_edge = np.array( time_edge ) + good_start
print( time_edge )
qrt_pds = get_t_qrc( FD, time_edge, Qrs, Qzs, qr_map, qz_map, mask=mask, path=data_dir, uid = uidstr )
plot_qrt_pds( qrt_pds, time_edge, qz_index = 0, uid = uidstr, path = data_dir )
In [72]:
#np.pi*2/0.022/10
In [73]:
if scat_geometry =='gi_saxs':
if run_profile_plot:
xcorners= [ 1100, 1250, 1250, 1100 ]
ycorners= [ 850, 850, 950, 950 ]
waterfall_roi_size = [ xcorners[1] - xcorners[0], ycorners[2] - ycorners[1] ]
waterfall_roi = create_rectangle_mask( avg_img, xcorners, ycorners )
#show_img( waterfall_roi * avg_img, aspect=1,vmin=.001, vmax=1, logs=True, )
wat = cal_waterfallc( FD, waterfall_roi, qindex= 1, bin_waterfall=True,
waterfall_roi_size = waterfall_roi_size,save =True, path=data_dir, uid=uidstr)
In [74]:
if scat_geometry =='gi_saxs':
if run_profile_plot:
plot_waterfallc( wat, qindex=1, aspect=None, vmin=1, vmax= np.max( wat), uid=uidstr, save =True,
path=data_dir, beg= FD.beg)
In [75]:
if scat_geometry =='gi_saxs':
show_qzr_roi( avg_img, roi_mask, inc_x0, ticks[:4], alpha=0.5, save=True, path=data_dir, uid=uidstr )
## Get 1D Curve (Q||-intensity¶)
qr_1d_pds = cal_1d_qr( avg_img, Qr, Qz, qr_map, qz_map, inc_x0= None, mask=mask, setup_pargs=setup_pargs )
plot_qr_1d_with_ROI( qr_1d_pds, qr_center=np.unique( np.array(list( qval_dict.values() ) )[:,0] ),
loglog=True, save=True, uid=uidstr, path = data_dir)
In [76]:
if scat_geometry =='gi_waxs':
#badpixel = np.where( avg_img[:600,:] >=300 )
#roi_mask[badpixel] = 0
show_ROI_on_image( avg_img, roi_mask, label_on = True, alpha=.5,
save=True, path=data_dir, uid=uidstr, vmin=0.1, vmax=500)
In [77]:
qind, pixelist = roi.extract_label_indices(roi_mask)
noqs = len(np.unique(qind))
print(noqs)
In [78]:
nopr = np.bincount(qind, minlength=(noqs+1))[1:]
nopr
Out[78]:
In [79]:
roi_inten = check_ROI_intensity( avg_img, roi_mask, ring_number= 1, uid =uidstr ) #roi starting from 1
In [80]:
qth_interest = 2 #the second ring. #qth_interest starting from 1
if scat_geometry =='saxs' or scat_geometry =='gi_waxs' or scat_geometry == 'gi_saxs':
if run_waterfall:
wat = cal_waterfallc( FD, roi_mask, qindex= qth_interest, save =True, path=data_dir, uid=uidstr)
plot_waterfallc( wat, qth_interest, aspect= None, vmin=1e-1, vmax= wat.max(), uid=uidstr, save =True,
path=data_dir, beg= FD.beg, cmap = cmap_vge )
In [81]:
#show_img(roi_mask)
In [82]:
ring_avg = None
if run_t_ROI_Inten:
times_roi, mean_int_sets = cal_each_ring_mean_intensityc(FD, roi_mask, timeperframe = None, multi_cor= True )#False )#True )
mean_int_setsF = np.zeros( [FD.end, mean_int_sets.shape[1] ] )
mean_int_setsF[FD.beg:FD.end] =mean_int_sets
plot_each_ring_mean_intensityc( times_roi, mean_int_sets, uid = uidstr, save=True, path=data_dir )
roi_avg = np.average( mean_int_sets, axis=0)
In [83]:
if run_get_mass_center:
cx, cy = get_mass_center_one_roi(FD, roi_mask, roi_ind=25)
In [84]:
if run_get_mass_center:
fig,ax=plt.subplots(2)
plot1D( cx, m='o', c='b',ax=ax[0], legend='mass center-refl_X',
ylim=[940, 960], ylabel='posX (pixel)')
plot1D( cy, m='s', c='r',ax=ax[1], legend='mass center-refl_Y',
ylim=[1540, 1544], xlabel='frames',ylabel='posY (pixel)')
Note : Enter the number of buffers for Muliti tau one time correlation number of buffers has to be even. More details in https://github.com/scikit-beam/scikit-beam/blob/master/skbeam/core/correlation.py
In [85]:
good_start_g2 = FD.beg
good_end_g2 = FD.end
define_good_series = False
#define_good_series = True
if define_good_series:
good_start_g2 = 1000 #399
good_end_g2 = 2000 #1486
FD = Multifile(filename, beg = good_start_g2, end = good_end_g2) #end=1000)
uid_ = uidstr + '_fra_%s_%s'%(FD.beg, FD.end)
print( uid_ )
In [86]:
#use_SG = False # True #False #True
#use_SG_bin_frames = True
In [87]:
if use_sqnorm:#for transmision SAXS
norm = get_pixelist_interp_iq( qp_saxs, iq_saxs, roi_mask, center)
print('Using circular average in the normalization of G2 for SAXS scattering.')
elif use_SG:#for Gi-SAXS or WAXS
avg_imgf = sgolay2d( avg_img, window_size= 11, order= 5) * mask
norm=np.ravel(avg_imgf)[pixelist]
print('Using smoothed image by SavitzkyGolay filter in the normalization of G2.')
elif use_SG_bin_frames:#for Gi-SAXS or WAXS
avg_imgf = sgolay2d( avg_img, window_size= 11, order= 5) * mask
norm_avg=np.ravel(avg_imgf)[pixelist]
#print('Using smoothed image by SavitzkyGolay filter in the normalization of G2.')
bins_number = 1
norm = get_SG_norm( FD, pixelist, bins=bins_number, mask=mask, window_size= 11, order=5 )
print('Using smoothed bined: %s frames by SavitzkyGolay filter in the normalization of G2.'%bins_number)
else:
norm= None
print('Using simple (average) normalization of G2.')
if use_imgsum_norm:
imgsum_ = imgsum
print('Using frame total intensity for intensity normalization in g2 calculation.')
else:
imgsum_ = None
import time
In [88]:
#np.save( data_dir + 'norm_SG_per_frame',norm)
#norm = np.load( data_dir + 'norm_SG_per_frame.npy' )
In [89]:
#%run -i ~/pyCHX_link/pyCHX/chx_correlationc.py
In [90]:
if run_one_time:
t0 = time.time()
if cal_g2_error:
g2,lag_steps,g2_err = cal_g2p(FD,roi_mask,bad_frame_list,good_start, num_buf = 8,
num_lev= None,imgsum= imgsum_, norm=norm, cal_error= True )
else:
g2,lag_steps = cal_g2p(FD,roi_mask,bad_frame_list,good_start, num_buf = 8,
num_lev= None,imgsum= imgsum_, norm=norm, cal_error= False )
run_time(t0)
In [91]:
lag_steps = lag_steps[:g2.shape[0]]
g2.shape[1]
Out[91]:
In [92]:
if run_one_time:
taus = lag_steps * timeperframe
try:
g2_pds = save_g2_general( g2, taus=taus,qr= np.array( list( qval_dict.values() ) )[:g2.shape[1],0],
qz = np.array( list( qval_dict.values() ) )[:g2.shape[1],1],
uid=uid_+'_g2.csv', path= data_dir, return_res=True )
except:
g2_pds = save_g2_general( g2, taus=taus,qr= np.array( list( qval_dict.values() ) )[:g2.shape[1],0],
uid=uid_+'_'+q_mask_name+'_g2.csv', path= data_dir, return_res=True )
if cal_g2_error:
try:
g2_err_pds = save_g2_general( g2_err, taus=taus,qr= np.array( list( qval_dict.values() ) )[:g2.shape[1],0],
qz = np.array( list( qval_dict.values() ) )[:g2.shape[1],1],
uid=uid_+'_g2_err.csv', path= data_dir, return_res=True )
except:
g2_err_pds = save_g2_general( g2_err, taus=taus,qr= np.array( list( qval_dict.values() ) )[:g2.shape[1],0],
uid=uid_+'_'+q_mask_name+'_g2_err.csv', path= data_dir, return_res=True )
In [ ]:
In [93]:
if run_one_time:
try:
g2_fit_result, taus_fit, g2_fit = get_g2_fit_general( g2, taus,
function = fit_g2_func, vlim=[0.95, 1.05], fit_range= None,
fit_variables={'baseline':False, 'beta': True, 'alpha':True,'relaxation_rate':True,},
guess_values={'baseline':1.0,'beta': 0.03,'alpha':1.0,'relaxation_rate':0.0005},
guess_limits = dict( baseline =[.9, 1.3], alpha=[0, 2],
beta = [0, 1], relaxation_rate= [1e-7, 1000]) ,)
g2_fit_paras = save_g2_fit_para_tocsv(g2_fit_result, filename= uid_ +'_'+q_mask_name +'_g2_fit_paras.csv', path=data_dir )
g2_fit_failed = False
except:
g2_fit_failed = True
In [94]:
scat_geometry_
Out[94]:
In [95]:
if run_one_time:
if cal_g2_error:
g2_fit_err = np.zeros_like(g2_fit)
if not g2_fit_failed:
plot_g2_general( g2_dict={1:g2, 2:g2_fit}, taus_dict={1:taus, 2:taus_fit},
vlim=[0.95, 1.05], g2_err_dict= {1:g2_err, 2: g2_fit_err},
qval_dict = dict(itertools.islice(qval_dict.items(),g2.shape[1])), fit_res= g2_fit_result, geometry= scat_geometry_,filename= uid_+'_g2',
path= data_dir, function= fit_g2_func, ylabel='g2', append_name= '_fit')
else:
plot_g2_general( g2_dict={1:g2, }, taus_dict={1:taus,},
vlim=[0.95, 1.05], g2_err_dict= {1:g2_err},
qval_dict = dict(itertools.islice(qval_dict.items(),g2.shape[1])),
fit_res= None, geometry= scat_geometry_,filename= uid_+'_g2',
path= data_dir, function= fit_g2_func, ylabel='g2', append_name= '')
else:
if not g2_fit_failed:
plot_g2_general( g2_dict={1:g2, 2:g2_fit}, taus_dict={1:taus, 2:taus_fit}, vlim=[0.95, 1.05],
qval_dict = dict(itertools.islice(qval_dict.items(),g2.shape[1])), fit_res= g2_fit_result, geometry= scat_geometry_,filename= uid_+'_g2',
path= data_dir, function= fit_g2_func, ylabel='g2', append_name= '_fit')
else:
plot_g2_general( g2_dict={1:g2}, taus_dict={1:taus}, vlim=[0.95, 1.05],
qval_dict = dict(itertools.islice(qval_dict.items(),g2.shape[1])), fit_res= None, geometry= scat_geometry_,filename= uid_+'_g2',
path= data_dir, function= fit_g2_func, ylabel='g2', append_name= '')
In [ ]:
In [ ]:
In [96]:
if run_one_time:
try:
if True:
fs, fe = 0, 8
#fs,fe=0, 6
qval_dict_ = {k:qval_dict[k] for k in list(qval_dict.keys())[fs:fe] }
D0, qrate_fit_res = get_q_rate_fit_general( qval_dict_, g2_fit_paras['relaxation_rate'][fs:fe],
geometry= scat_geometry_ )
plot_q_rate_fit_general( qval_dict_, g2_fit_paras['relaxation_rate'][fs:fe], qrate_fit_res,
geometry= scat_geometry_,uid=uid_ , path= data_dir )
elif False:
D0, qrate_fit_res = get_q_rate_fit_general( qval_dict, g2_fit_paras['relaxation_rate'],
fit_range=[0, 26], geometry= scat_geometry_ )
plot_q_rate_fit_general( qval_dict, g2_fit_paras['relaxation_rate'], qrate_fit_res,
geometry= scat_geometry_,uid=uid_ ,
show_fit=False, path= data_dir, plot_all_range=False)
else:
##plot log_q, log_rate
plot_q_rate_general( qval_dict, g2_fit_paras['relaxation_rate'],
geometry= scat_geometry_,uid=uid_ , path= data_dir )
plot_q_g2fitpara_general( qval_dict, g2_fit_paras,
geometry= scat_geometry_,uid=uid_ , path= data_dir )
except:
pass
In [97]:
#plot1D( x= qr, y=g2_fit_paras['beta'], ls='-', m = 'o', c='b', ylabel=r'$\beta$', xlabel=r'$Q( \AA^{-1} ) $' )
In [98]:
define_good_series = False
#define_good_series = True
if define_good_series:
good_start = 0
good_end = 9000
FD = Multifile(filename, beg = good_start, end = good_end)
uid_ = uidstr + '_fra_%s_%s'%(FD.beg, FD.end)
print( uid_ )
In [99]:
data_pixel = None
if run_two_time:
data_pixel = Get_Pixel_Arrayc( FD, pixelist, norm= norm, ).get_data()
### With another normalizaiton by using roi_intensity for each frame
#data_pixel = Get_Pixel_Arrayc( FD, pixelist, norm= norm, mean_int_sets = mean_int_setsF,qind=qind ).get_data()
In [100]:
import time
t0=time.time()
g12b=None
if run_two_time:
g12b = auto_two_Arrayc( data_pixel, roi_mask, index = None )
if run_dose:
np.save( data_dir + 'uid=%s_g12b'%uid, g12b)
run_time( t0 )
In [101]:
if run_two_time:
show_C12( np.nan_to_num( g12b ), q_ind= 2, qlabel=dict(itertools.islice(qval_dict.items(),g12b.shape[2])),
#N1= FD.beg,logs=False, N2=min( FD.end,10000),
vmin= 1.00, vmax=1.1, timeperframe=timeperframe,save=True,
path= data_dir, uid = uid_ ,cmap=plt.cm.jet)#cmap=cmap_albula)
In [102]:
multi_tau_steps = True
if run_two_time:
if lag_steps is None:
num_bufs=8
noframes = FD.end - FD.beg
num_levels = int(np.log( noframes/(num_bufs-1))/np.log(2) +1) +1
tot_channels, lag_steps, dict_lag = multi_tau_lags(num_levels, num_bufs)
max_taus= lag_steps.max()
#max_taus= lag_steps.max()
max_taus = Nimg
t0=time.time()
#tausb = np.arange( g2b.shape[0])[:max_taus] *timeperframe
if multi_tau_steps:
lag_steps_ = lag_steps[ lag_steps <= g12b.shape[0] ]
g2b = get_one_time_from_two_time(g12b)[lag_steps_]
tausb = lag_steps_ *timeperframe
else:
tausb = (np.arange( g12b.shape[0]) *timeperframe)[:-200]
g2b = (get_one_time_from_two_time(g12b))[:-200]
run_time(t0)
g2b_pds = save_g2_general( g2b, taus=tausb, qr= np.array( list( qval_dict.values() ) )[:g2.shape[1],0],
qz=None, uid=uid_+'_'+q_mask_name+'_g2b.csv', path= data_dir, return_res=True )
In [103]:
if run_two_time:
g2b_fit_result, tausb_fit, g2b_fit = get_g2_fit_general( g2b, tausb,
function = fit_g2_func, vlim=[0.95, 1.05], fit_range= None,
fit_variables={'baseline':False, 'beta': True, 'alpha':True,'relaxation_rate':True},
guess_values={'baseline':1.0,'beta': 0.15,'alpha':1.0,'relaxation_rate':1e-3,},
guess_limits = dict( baseline =[1, 1.8], alpha=[0, 2],
beta = [0, 1], relaxation_rate= [1e-8, 5000]) )
g2b_fit_paras = save_g2_fit_para_tocsv(g2b_fit_result, filename= uid_ +'_'+q_mask_name+'_g2b_fit_paras.csv', path=data_dir )
In [104]:
#plot1D( x = tausb[1:], y =g2b[1:,0], ylim=[0.95, 1.46], xlim = [0.0001, 10], m='', c='r', ls = '-',
# logx=True, title='one_time_corelation', xlabel = r"$\tau $ $(s)$", )
In [105]:
if run_two_time:
plot_g2_general( g2_dict={1:g2b, 2:g2b_fit}, taus_dict={1:tausb, 2:tausb_fit}, vlim=[0.95, 1.05],
qval_dict=dict(itertools.islice(qval_dict.items(),g2.shape[1])), fit_res= g2b_fit_result, geometry=scat_geometry_,filename=uid_+'_g2',
path= data_dir, function= fit_g2_func, ylabel='g2', append_name= '_b_fit')
In [106]:
if run_two_time:
D0b, qrate_fit_resb = get_q_rate_fit_general( dict(itertools.islice(qval_dict.items(),g2.shape[1])), g2b_fit_paras['relaxation_rate'],
fit_range=[0, 10], geometry= scat_geometry_ )
In [107]:
#qval_dict, g2b_fit_paras['relaxation_rate']
In [108]:
#%run ~/pyCHX_link/pyCHX/chx_generic_functions.py
In [109]:
if run_two_time:
try:
if True:
fs, fe = 0, 8
#fs,fe=0, 6
qval_dict_ = {k:qval_dict[k] for k in list(qval_dict.keys())[fs:fe] }
D0, qrate_fit_res = get_q_rate_fit_general( qval_dict_, g2b_fit_paras['relaxation_rate'][fs:fe],
geometry= scat_geometry_ )
plot_q_rate_fit_general( qval_dict_, g2b_fit_paras['relaxation_rate'][fs:fe], qrate_fit_res,
geometry= scat_geometry_,uid=uid_ +'_two_time' , path= data_dir )
elif False:
D0, qrate_fit_res = get_q_rate_fit_general( qval_dict, g2b_fit_paras['relaxation_rate'],
fit_range=[0, 26], geometry= scat_geometry_ )
plot_q_rate_fit_general( qval_dict, g2_fit_paras['relaxation_rate'], qrate_fit_res,
geometry= scat_geometry_,uid=uid_ +'_two_time' ,
show_fit=False, path= data_dir, plot_all_range=False)
else:
##plot log_q, log_rate
plot_q_rate_general( qval_dict, g2b_fit_paras['relaxation_rate'],
geometry= scat_geometry_,uid=uid_ +'_two_time' , path= data_dir )
plot_q_g2fitpara_general( qval_dict, g2b_fit_paras,
geometry= scat_geometry_,uid=uid_ +'_two_time' , path= data_dir )
except:
pass
In [110]:
if run_two_time and run_one_time:
plot_g2_general( g2_dict={1:g2, 2:g2b}, taus_dict={1:taus, 2:tausb},vlim=[0.99, 1.007],
qval_dict=dict(itertools.islice(qval_dict.items(),g2.shape[1])), g2_labels=['from_one_time', 'from_two_time'],
geometry=scat_geometry_,filename=uid_+'_g2_two_g2', path= data_dir, ylabel='g2', )
In [111]:
#run_dose = True
In [112]:
if run_dose:
get_two_time_mulit_uids( [uid], roi_mask, norm= norm, bin_frame_number=1,
path= data_dir0, force_generate=False, compress_path = cmp_path + '/' )
In [113]:
try:
print( md['transmission'] )
except:
md['transmission'] =1
In [114]:
exposuretime
Out[114]:
In [115]:
if run_dose:
N = len(imgs)
print(N)
#exposure_dose = md['transmission'] * exposuretime* np.int_([ N/16, N/8, N/4 ,N/2, 3*N/4, N*0.99 ])
exposure_dose = md['transmission'] * exposuretime* np.int_([ N/8, N/4 ,N/2, 3*N/4, N*0.99 ])
print( exposure_dose )
In [116]:
if run_dose:
taus_uids, g2_uids = get_series_one_time_mulit_uids( [ uid ], qval_dict, good_start=good_start,
path= data_dir0, exposure_dose = exposure_dose, num_bufs =8, save_g2= False,
dead_time = 0, trans = [ md['transmission'] ] )
In [117]:
if run_dose:
plot_dose_g2( taus_uids, g2_uids, ylim=[1.0, 1.2], vshift= 0.00,
qval_dict = qval_dict, fit_res= None, geometry= scat_geometry_,
filename= '%s_dose_analysis'%uid_,
path= data_dir, function= None, ylabel='g2_Dose', g2_labels= None, append_name= '' )
In [118]:
if run_dose:
qth_interest = 5
plot_dose_g2( taus_uids, g2_uids, qth_interest= qth_interest, ylim=[0.98, 1.4], vshift= 0.00,
qval_dict = qval_dict, fit_res= None, geometry= scat_geometry_,
filename= '%s_dose_analysis'%uidstr,
path= data_dir, function= None, ylabel='g2_Dose', g2_labels= None, append_name= '' )
In [119]:
if run_four_time:
t0=time.time()
g4 = get_four_time_from_two_time(g12b, g2=g2b)[:int(max_taus)]
run_time(t0)
In [120]:
if run_four_time:
taus4 = np.arange( g4.shape[0])*timeperframe
g4_pds = save_g2_general( g4, taus=taus4, qr=np.array( list( qval_dict.values() ) )[:,0],
qz=None, uid=uid_ +'_g4.csv', path= data_dir, return_res=True )
In [121]:
if run_four_time:
plot_g2_general( g2_dict={1:g4}, taus_dict={1:taus4},vlim=[0.95, 1.05], qval_dict=qval_dict, fit_res= None,
geometry=scat_geometry_,filename=uid_+'_g4',path= data_dir, ylabel='g4')
In [122]:
if run_xsvs:
max_cts = get_max_countc(FD, roi_mask )
#max_cts = 15 #for eiger 500 K
qind, pixelist = roi.extract_label_indices( roi_mask )
noqs = len( np.unique(qind) )
nopr = np.bincount(qind, minlength=(noqs+1))[1:]
#time_steps = np.array( utils.geometric_series(2, len(imgs) ) )
time_steps = [0,1] #only run the first two levels
num_times = len(time_steps)
times_xsvs = exposuretime + (2**( np.arange( len(time_steps) ) ) -1 ) * timeperframe
print( 'The max counts are: %s'%max_cts )
In [123]:
if run_xsvs:
if roi_avg is None:
times_roi, mean_int_sets = cal_each_ring_mean_intensityc(FD, roi_mask, timeperframe = None, )
roi_avg = np.average( mean_int_sets, axis=0)
t0=time.time()
spec_bins, spec_his, spec_std, spec_sum = xsvsp( FD, np.int_(roi_mask), norm=None,
max_cts=int(max_cts+2), bad_images=bad_frame_list, only_two_levels=True )
spec_kmean = np.array( [roi_avg * 2**j for j in range( spec_his.shape[0] )] )
run_time(t0)
spec_pds = save_bin_his_std( spec_bins, spec_his, spec_std, filename=uid_+'_spec_res.csv', path=data_dir )
In [124]:
if run_xsvs:
ML_val, KL_val,K_ = get_xsvs_fit( spec_his, spec_sum, spec_kmean,
spec_std, max_bins=2, fit_range=[1,60], varyK= False )
#print( 'The observed average photon counts are: %s'%np.round(K_mean,4))
#print( 'The fitted average photon counts are: %s'%np.round(K_,4))
print( 'The difference sum of average photon counts between fit and data are: %s'%np.round(
abs(np.sum( spec_kmean[0,:] - K_ )),4))
print( '#'*30)
qth= 0
print( 'The fitted M for Qth= %s are: %s'%(qth, ML_val[qth]) )
print( K_[qth])
print( '#'*30)
In [ ]:
In [125]:
if run_xsvs:
qr = [qval_dict[k][0] for k in list(qval_dict.keys()) ]
plot_xsvs_fit( spec_his, ML_val, KL_val, K_mean = spec_kmean, spec_std=spec_std,
xlim = [0,10], vlim =[.9, 1.1],
uid=uid_, qth= qth_interest, logy= True, times= times_xsvs, q_ring_center=qr, path=data_dir)
plot_xsvs_fit( spec_his, ML_val, KL_val, K_mean = spec_kmean, spec_std = spec_std,
xlim = [0,15], vlim =[.9, 1.1],
uid=uid_, qth= None, logy= True, times= times_xsvs, q_ring_center=qr, path=data_dir )
In [126]:
if run_xsvs:
contrast_factorL = get_contrast( ML_val)
spec_km_pds = save_KM( spec_kmean, KL_val, ML_val, qs=qr, level_time=times_xsvs, uid=uid_, path = data_dir )
#spec_km_pds
In [127]:
if run_xsvs:
plot_g2_contrast( contrast_factorL, g2b, times_xsvs, tausb, qr,
vlim=[0.8,1.2], qth = qth_interest, uid=uid_,path = data_dir, legend_size=14)
plot_g2_contrast( contrast_factorL, g2b, times_xsvs, tausb, qr,
vlim=[0.8,1.2], qth = None, uid=uid_,path = data_dir, legend_size=4)
In [128]:
#from chxanalys.chx_libs import cmap_vge, cmap_albula, Javascript
In [129]:
md['mask_file']= mask_path + mask_name
md['roi_mask_file']= fp
md['mask'] = mask
#md['NOTEBOOK_FULL_PATH'] = data_dir + get_current_pipeline_fullpath(NFP).split('/')[-1]
md['good_start'] = good_start
md['bad_frame_list'] = bad_frame_list
md['avg_img'] = avg_img
md['roi_mask'] = roi_mask
md['setup_pargs'] = setup_pargs
if scat_geometry == 'gi_saxs':
md['Qr'] = Qr
md['Qz'] = Qz
md['qval_dict'] = qval_dict
md['beam_center_x'] = inc_x0
md['beam_center_y']= inc_y0
md['beam_refl_center_x'] = refl_x0
md['beam_refl_center_y'] = refl_y0
elif scat_geometry == 'gi_waxs':
md['beam_center_x'] = center[1]
md['beam_center_y']= center[0]
else:
md['qr']= qr
#md['qr_edge'] = qr_edge
md['qval_dict'] = qval_dict
md['beam_center_x'] = center[1]
md['beam_center_y']= center[0]
md['beg'] = FD.beg
md['end'] = FD.end
md['qth_interest'] = qth_interest
md['metadata_file'] = data_dir + 'uid=%s_md.pkl'%uid
psave_obj( md, data_dir + 'uid=%s_md.pkl'%uid ) #save the setup parameters
save_dict_csv( md, data_dir + 'uid=%s_md.csv'%uid, 'w')
###############Only work for Randy
try:
md['beg_OneTime'] = good_start_g2
md['end_OneTime'] = good_end_g2
print(md['beg_OneTime'])
print(md['end_OneTime'])
except:
pass
################################
Exdt = {}
if scat_geometry == 'gi_saxs':
for k,v in zip( ['md', 'roi_mask','qval_dict','avg_img','mask','pixel_mask', 'imgsum', 'bad_frame_list', 'qr_1d_pds'],
[md, roi_mask, qval_dict, avg_img,mask,pixel_mask, imgsum, bad_frame_list, qr_1d_pds] ):
Exdt[ k ] = v
elif scat_geometry == 'saxs':
for k,v in zip( ['md', 'q_saxs', 'iq_saxs','iqst','qt','roi_mask','qval_dict','avg_img','mask','pixel_mask', 'imgsum', 'bad_frame_list'],
[md, q_saxs, iq_saxs, iqst, qt,roi_mask, qval_dict, avg_img,mask,pixel_mask, imgsum, bad_frame_list] ):
Exdt[ k ] = v
elif scat_geometry == 'gi_waxs':
for k,v in zip( ['md', 'roi_mask','qval_dict','avg_img','mask','pixel_mask', 'imgsum', 'bad_frame_list'],
[md, roi_mask, qval_dict, avg_img,mask,pixel_mask, imgsum, bad_frame_list] ):
Exdt[ k ] = v
if run_waterfall:Exdt['wat'] = wat
if run_t_ROI_Inten:Exdt['times_roi'] = times_roi;Exdt['mean_int_sets']=mean_int_sets
if run_one_time:
if run_invariant_analysis:
for k,v in zip( ['taus','g2','g2_fit_paras', 'time_stamp','invariant'], [taus,g2,g2_fit_paras,time_stamp,invariant] ):Exdt[ k ] = v
else:
for k,v in zip( ['taus','g2','g2_fit_paras' ], [taus,g2,g2_fit_paras ] ):Exdt[ k ] = v
if run_two_time:
for k,v in zip( ['tausb','g2b','g2b_fit_paras', 'g12b'], [tausb,g2b,g2b_fit_paras,g12b] ):Exdt[ k ] = v
#for k,v in zip( ['tausb','g2b','g2b_fit_paras', ], [tausb,g2b,g2b_fit_paras] ):Exdt[ k ] = v
if run_dose:
for k,v in zip( [ 'taus_uids', 'g2_uids' ], [taus_uids, g2_uids] ):Exdt[ k ] = v
if run_four_time:
for k,v in zip( ['taus4','g4'], [taus4,g4] ):Exdt[ k ] = v
if run_xsvs:
for k,v in zip( ['spec_kmean','spec_pds','times_xsvs','spec_km_pds','contrast_factorL'],
[ spec_kmean,spec_pds,times_xsvs,spec_km_pds,contrast_factorL] ):Exdt[ k ] = v
In [130]:
export_xpcs_results_to_h5( 'uid=%s_%s_Res.h5'%(md['uid'],q_mask_name), data_dir, export_dict = Exdt )
#extract_dict = extract_xpcs_results_from_h5( filename = 'uid=%s_Res.h5'%md['uid'], import_dir = data_dir )
In [131]:
from pyCHX.chx_generic_functions import append_txtfile
outp = '/nsls2/xf11id1/analysis/2019_3/%s/'%username
append_txtfile( filename = outp + username + '.lst',
data = np.array( [ outp + 'Results/' + md['uid'] + '/uid=%s__Res.h5'%md['uid']]))
In [132]:
#g2npy_filename = data_dir + '/' + 'uid=%s_g12b.npy'%uid
#print(g2npy_filename)
#if os.path.exists( g2npy_filename):
# print('Will delete this file=%s.'%g2npy_filename)
# os.remove( g2npy_filename )
In [133]:
#extract_dict = extract_xpcs_results_from_h5( filename = 'uid=%s_Res.h5'%md['uid'], import_dir = data_dir )
In [134]:
#extract_dict = extract_xpcs_results_from_h5( filename = 'uid=%s_Res.h5'%md['uid'], import_dir = data_dir )
In [135]:
#fp = data_dir + 'uid=4c747cec-3708-404b-a5d7-b909e9659b83__Res.h5'
#f = h5py.File(fp , 'r')
In [ ]:
In [ ]:
In [136]:
pdf_out_dir = os.path.join('/XF11ID/analysis/', CYCLE, username, 'Results/')
pdf_filename = "XPCS_Analysis_Report2_for_uid=%s%s%s.pdf"%(uid,pdf_version,q_mask_name)
if run_xsvs:
pdf_filename = "XPCS_XSVS_Analysis_Report_for_uid=%s%s%s.pdf"%(uid,pdf_version,q_mask_name)
In [137]:
#%run /home/yuzhang/pyCHX_link/pyCHX/Create_Report.py
In [138]:
make_pdf_report( data_dir, uid, pdf_out_dir, pdf_filename, username,
run_fit_form,run_one_time, run_two_time, run_four_time, run_xsvs, run_dose,
report_type= scat_geometry, report_invariant= run_invariant_analysis,
md = md )
In [6]:
update_olog_uid_with_file?
In [140]:
if att_pdf_report:
lops=0
while lops<5:
try:
os.environ['HTTPS_PROXY'] = 'https://proxy:8888'
os.environ['no_proxy'] = 'cs.nsls2.local,localhost,127.0.0.1'
update_olog_uid_with_file( uid[:6], text='Add XPCS Analysis PDF Report',
filename=pdf_out_dir + pdf_filename, append_name='_R1' )
lops=5
except:
print("Failed to attach PDF report to Olog...try again in 10s")
time.sleep(10)
lops+=1
if lops == 5:
print("Failed final attempt to attach PDF to Olog!")
In [141]:
#save_oavs= False #True
In [142]:
if save_oavs:
try:
os.environ['HTTPS_PROXY'] = 'https://proxy:8888'
os.environ['no_proxy'] = 'cs.nsls2.local,localhost,127.0.0.1'
update_olog_uid_with_file( uid[:6], text='Add OVA images',
filename= data_dir + 'uid=%s_OVA_images.png'%uid, append_name='_img' )
except:
pass
In [ ]:
In [ ]: