NuBS准弹散射谱的数据规约

博士后出站报告附录

# -*- coding: utf-8 -*-
"""
Created on Tue Dec 21 10:50:15 2021

@author: zhouhuibin
"""

import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
import numpy as np
import math
import os
from zhbcode import selffunction as sf
from zhbcode.find_data import load_McStasdata,load_weidata
from tkinter import _flatten # 将二维列表转换为一维列表
from scipy import interpolate

# Physics constant
m_n = 1.6749286e-27 # unit kg, neutron mass
h_P = 6.62607015e-34 # unit J s, planck constant
e = 1.602176634e-19

# incident beam parameters
L_i = L_1 = 90 # unit m   

# Analyzer parameter
d = 3.13545e-10 # unit m, Spaing of Si-111 plane
lambda_off = 2*d # wavelength when theta is 90 deg;

t_0 = 46.887e-6 #unit s, time delay for 6.267 AA 

h_d = 0.15 # unit m, the vertical offset of the detector center;
R = 2.5 # unit m, radius of analyzer sphere
y_off = np.linspace(-0.2,0.2,41)
h_0 = y_off+h_d
tof = np.linspace(0.14900, 0.150800,200)

maindir = r'D:\001_CSNS\data\NBSS_data_reduction\TOF_DATA\Python\20230202_Spectrum_Sample_33_Analyser_1mm'
datapath = sf.mkdir(maindir+r'\Norm Data Output')
figpath = sf.mkdir(maindir+r'\figs')

rawdatapath = maindir+r'\raw_data'
namelist = os.listdir(rawdatapath)

data_SecondGeo_file = maindir+r'\Second_Spec_data\223_alldata_0.00cm.dat'
data_SecondGeo = load_weidata(data_SecondGeo_file,'0.00cm')
theta_B_raw = data_SecondGeo[:,2]
y_off_raw = data_SecondGeo[:,6]
L_f_raw = data_SecondGeo[:,9]
L_2_raw = data_SecondGeo[:,10]
phi_deg_raw = data_SecondGeo[:,3]

# data_t_0_file = maindir+r'\t0_data\Energy-t_select.dat'
# data_t_0 = load_weidata(data_t_0_file,'DPHM-BL10')
# Energy_i_raw = data_t_0[:,0]
# t_0_raw = data_t_0[:,1]
# t_0_min = min(t_0_raw)
# t_0_max = max(t_0_raw)

theta_B_degs = interpolate.interp1d(y_off_raw, theta_B_raw, kind='linear')(y_off)
theta_B_rads = np.radians(theta_B_degs)
L_fs = interpolate.interp1d(y_off_raw, L_f_raw, kind='linear')(y_off)
L_2s = interpolate.interp1d(y_off_raw, L_2_raw, kind='linear')(y_off)
phis_deg = interpolate.interp1d(y_off_raw, phi_deg_raw, kind='linear')(y_off)
phis_rad = np.radians(phis_deg)
phi_fs_deg = np.linspace(0,160,33)
phi_fs_rad = np.radians(phi_fs_deg)

Q_all = []
E_all = []
I_all = []
posi_pixel = y_off

name_para = []
unit_para = []
comment_para = []
data_para = []

name_E0 = []
unit_E0 = []
comment_E0 = []
data_E0 = []

name_FWHM = []
unit_FWHM = []
comment_FWHM = []
data_FWHM = []

name_tcenter = ['h_0',]
unit_tcenter = ['m',]
comment_tcenter = ['none',]
data_tcenter = [h_0]

for i0 in (range(len(namelist))):
    data = load_McStasdata(os.path.join(rawdatapath,namelist[i0]), '#')
    print (namelist[i0],data.shape)
    E_is = [] # 存放每一个探测器的所有initial Energy
    E_fs = []   # 存放每一个探测器的final Energy
    Es = [] #存放一个探测器的Energy transfer
    Is = [] #存放一个探测器的Intensity

    name_pixel = []
    unit_pixel = []
    comment_pixel = []
    data_pixel = []

    phi_f_rad = phi_fs_rad[i0]
    Q_eachDet = []
    E_0_eachDet = []
    fwhm_eachDet = []
    t_centers = []

    for i1 in range(data.shape[0]): #探测器TOF谱的行数，即像素点的个数
        cos_2theta = (R/L_2s[i1])*np.sin(phis_rad[i1])*np.cos(phi_f_rad)
        Q = (4*np.pi/6.267)*np.sin(math.acos(cos_2theta)/2) # unit 1/A
        Q_all.append(Q) #所有探测器所有像素点的Q
        t_single= []
        E_is_single = [] #存放每一个像素点的E_i
        E_fs_single = [] #存放每一个像素点的E_f
        Es_single = []  # 存放每一个像素点的Energy transfer
        Is_single = []  # 存放每一个像素点的Intensity
        for i2 in range(data.shape[1]): # 探测器TOF谱的列数，即飞行时间点的个数
            I = data[i1,i2] # 第i1个像素，第i2个飞行时间点的Intensity
            t = tof[i2] # 第i2个飞行时间点对应的飞行时间
            theta_B_rad = theta_B_rads[i1] # 第i1个像素点对应的theta_B
            L_f = L_fs[i1] # 第i1个像素点对应的飞行距离
            lambda_f = 2*d*np.sin(theta_B_rad)
            v_f = h_P/(lambda_f*m_n)
            t_2 = L_f/v_f
            '''
            def Delta_E_i(t_0_inter): # unit of t_0_inter is us
                E_i_source = interpolate.interp1d(t_0_raw, Energy_i_raw, kind='linear')(t_0_inter) # unit is meV
                E_i_tof = 0.5*m_n*(L_1/(t-t_2-t_0_inter*1e-6))**2/(e*1e-3) # unit is meV
                delta_E_i = E_i_source - E_i_tof
                return delta_E_i

            t_0 = sf.biSection(t_0_min,t_0_max,1e-4,Delta_E_i)*1e-6 #unit is second
            '''
            E_i = 0.5*m_n*(L_1/(t-t_2-t_0))**2
            E_f = 0.5/m_n*(h_P/lambda_f)**2
            E = E_i - E_f

            E_is.append(E_i)
            E_fs.append(E_f)
            Es.append(E)
            Is.append(I)

            t_single.append(t)
            E_is_single.append(E_i)
            E_fs_single.append(E_f)
            Es_single.append(E)
            Is_single.append(I)
        t_single = np.array(t_single)
        E_is_single = np.array(E_is_single)
        E_fs_single = np.array(E_fs_single)
        Es_single = np.array(Es_single)
        Is_single = np.array(Is_single)

        E_break = -10*(1e-6*e)  # 切断多余的能量（可能有其它的峰）
        for i3 in range(len(Es_single)):
            if Es_single[i3] < E_break:
                index_break = i3
                break
        t_single= t_single[0:index_break+1]
        Es_single = Es_single[0:index_break+1]
        Is_single = Is_single[0:index_break+1]

        Es_grid = np.linspace(-10*(1e-6*e),25*(1e-6*e),201)

        # print (min(Es_single)/(1e-6*e),max(Es_single)/(1e-6*e))


        t_grid = interpolate.interp1d(Es_single, t_single, kind='linear')(Es_grid)
        Is_grid = interpolate.interp1d(Es_single, Is_single, kind='linear')(Es_grid)
        t_single = t_grid
        Es_single = Es_grid
        Is_single = Is_grid

        if max(Is_single) == 0:
            t_center = 0
        else:
            t_center = t_single[np.where(Is_single == max(Is_single))[0][0]]
        t_centers.append(t_center)

        # print(len(Es_single),len(Is_single))

        E_0_single,fwhm_single = sf.FWHM(Es_single,Is_single)
        E_0_eachDet.append(E_0_single)
        fwhm_eachDet.append(fwhm_single)
        Q_eachDet.append(Q)


        E_all.append(Es_single) #所有探测器的所有像素点的E
        I_all.append(Is_single) #所有探测器的所有像素点的I

        # TOF-pixel output
        name_pixel = name_pixel+['tof','Energy_transfer','Intensity',]
        unit_pixel = unit_pixel+['s','ueV','none',]
        comment_pixel = comment_pixel+[(f'{i1}-pixel-{round(y_off[i1],2)}m-{round(Q,2)}A-1',)*3]
        data_pixel = data_pixel+[t_single,Es_single/(1e-6*e),Is_single]

        # if i1 == 23: # plot第23个像素点的散射谱
        #     plt.cla()
        #     plt.clf()
        #     plt.plot(Es_single/(1e-6*e),Is_single,'ro')
        #     plt.xlabel('E(ueV)')
        #     plt.ylabel('I')
        #     plt.savefig(os.path.join(figpath,f'Is-single_vs_Es-single_{i0}-{i1}.png'))
        #     plt.show()

    # spectrum for each pixel     
    comment_pixel = _flatten(comment_pixel)
    data_pixel = np.transpose(data_pixel)
    filename = str(i0).rjust(3,'0')+f'_Spectrum_{i0}-Pixel_seprate.dat'
    sf.dataoutput(name_pixel,unit_pixel,comment_pixel,data_pixel,datapath+r'\Spectrum_pixel',filename)

    E_is = np.array(E_is)
    E_fs = np.array(E_fs)
    Es = np.array(Es)
    Is = np.array(Is)

    E_0_eachDet = np.array(E_0_eachDet)
    fwhm_eachDet = np.array(fwhm_eachDet)

    # spectrum for all  pixel for one detector
    name = ['E_i','E_f','Energy_transfer','Intensity']
    unit = ['ueV','ueV','ueV','none']
    comment = ['6.267AA','6.267AA','6.267AA','6.267AA']
    data = np.c_[E_is/(1e-6*e),E_fs/(1e-6*e),Es/(1e-6*e),Is]
    filename = str(i0).rjust(3,'0')+f'_Spectrum_{i0}.dat'
    sf.dataoutput(name,unit,comment,data,datapath+r'\Spectrum_pixel',filename)

    name_para = name_para+['Position','Q','E0','FWHM',]
    unit_para = unit_para+['m','A-1','ueV','ueV',]
    comment_para = comment_para+[(f'{phi_fs_deg[i0]}deg',)*4]
    data_para = data_para+[posi_pixel, Q_eachDet, E_0_eachDet/(1e-6*e),fwhm_eachDet/(1e-6*e)]

    name_E0 = name_E0+['E0',]
    unit_E0 = unit_E0+['ueV',]
    comment_E0 = comment_E0+[(f'{phi_fs_deg[i0]}deg',)*1]
    data_E0 = data_E0+[E_0_eachDet/(1e-6*e)]

    name_FWHM = name_FWHM+['FWHM',]
    unit_FWHM = unit_FWHM+['ueV',]
    comment_FWHM = comment_FWHM+[(f'{phi_fs_deg[i0]}deg',)*1]
    data_FWHM = data_FWHM+[fwhm_eachDet/(1e-6*e)]

    name_tcenter = name_tcenter + ['t_center',]
    unit_tcenter = unit_tcenter + ['s',]
    comment_tcenter = comment_tcenter + [f'{phi_fs_deg[i0]}deg',]
    data_tcenter = data_tcenter + [t_centers]

comment_para = _flatten(comment_para)
data_para = np.transpose(data_para)
filename = 'Parameter.dat'
sf.dataoutput(name_para,unit_para,comment_para,data_para,datapath,filename)

comment_E0 = _flatten(comment_E0)
data_E0 = np.transpose(data_E0)
filename = 'Parameter_E0.dat'
sf.dataoutput(name_E0,unit_E0,comment_E0,data_E0,datapath,filename)

comment_FWHM = _flatten(comment_FWHM)
data_FWHM = np.transpose(data_FWHM)
filename = 'Parameter_FWHM.dat'
sf.dataoutput(name_FWHM,unit_FWHM,comment_FWHM,data_FWHM,datapath,filename)

comment_tcenter = _flatten(comment_tcenter)
data_tcenter = np.transpose(data_tcenter)
filename = 't_center.dat'
sf.dataoutput(name_tcenter,unit_tcenter,comment_tcenter,data_tcenter,datapath,filename)

E_grid = np.linspace(0.9*min(Es_single),0.9*max(Es_single),200) #min(Es_single)为负值
Q_value = [0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8,2]
Q_space = 0.05
E_0s = []
fwhms = []

for i3 in range(len(Q_value)): # 合成每一个Q值的Spectrum
    Q_min = Q_value[i3]-Q_space
    Q_max = Q_value[i3]+Q_space
    I_grids = []

    name_Q = ['E']
    unit_Q = ['ueV']
    comment_Q = ['none']
    data_Q = [E_grid/(1e-6*e)]

    I_sum = np.zeros(len(E_grid))
    for i4 in range(len(Q_all)): 
        if Q_all[i4] >= Q_min and Q_all[i4] <= Q_max: #判断每一个像素的Q值是否在选定的Q值范围内
            I_grid = interpolate.interp1d(E_all[i4], I_all[i4], kind='linear')(E_grid) # 插值到同一个范围
            I_sum = np.array(I_sum)+np.array(I_grid)

            name_Q = name_Q+['I',]
            unit_Q = unit_Q+['none',]
            comment_Q = comment_Q+[f'{round(Q_all[i4],2)}A-1']
            data_Q = data_Q+[I_grid]

    name_Q = name_Q+['I_sum',]
    unit_Q = unit_Q+['none',]
    comment_Q = comment_Q+[f'{round(Q_value[i3],2)}A-1']
    data_Q = data_Q+[I_sum]

    comment_Q = _flatten(comment_Q)
    data_Q = np.transpose(data_Q)
    filename = str(i3).rjust(3,'0')+f'_{round(Q_value[i3],2)}A-1_Spectrum-Q.dat'
    sf.dataoutput(name_Q,unit_Q,comment_Q,data_Q,datapath+r'\Spectrum_Q',filename)

    E_0,fwhm = sf.FWHM(E_grid,I_sum)
    E_0s.append(E_0)
    fwhms.append(fwhm)
    print(i3)
E_0s = np.array(E_0s)
fwhms = np.array(fwhms)

name = _flatten(['Q_value','E_0','FWHM'])
unit = _flatten(['1/A','ueV','ueV'])
comment = _flatten([('Vanadium',)*3])
data = np.c_[Q_value,E_0s/(1e-6*e),fwhms/(1e-6*e)]
filename = str(i3+1).rjust(3,'0')+'_E_0_FWHM.dat'
sf.dataoutput(name,unit,comment,data,datapath+r'\Spectrum_Q',filename)