计算不同探测器管不同像素点处的动量转移大小

博士后出站报告附录

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

@author: zhouhuibin
"""

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

maindir = r'D:\CSNS\Summary\Data Treatment\python\Momentum resolution\005_Q-and-Q-Resolution'
datapath = sf.mkdir(maindir+r'\Norm Data Output')
figpath = sf.mkdir(maindir+r'\figs')

# Physics constant
m_n = 1.675e-27 # unit kg, neutron mass
h_P = 6.626e-34 # unit J s, planck constant

# incident beam parameters
L_i = L_1 = 90 # unit m
lambda_i_array = np.array([2,3,4,5,6.267,7,8,9,10])*1e-10 #unit m
theta_i_deg = 90 # unit deg
theta_i_rad = np.radians(theta_i_deg)
phi_i_deg = 0 # unit deg
phi_i_rad = np.radians(phi_i_deg)

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

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

theta_Dm_deg = np.linspace(-10,-160,16)
theta_Dp_deg = np.linspace(10,160,16)
theta_D_deg = np.concatenate([theta_Dm_deg,theta_Dp_deg]) # horizontal cover angle of the detectors

# Detector parameters
h_0_min = 0.129
h_0_max = 0.277
h_0s = np.linspace(h_0_min,h_0_max,round(100*(h_0_max-h_0_min))+1)
h_0_inters = []
for i0 in range(len(h_0s)-1): 
    h_0_inter = (h_0s[i0+1]+h_0s[i0])/2
    h_0_inters.append(h_0_inter)
# Parameters of secondary spectrometer
h_sample = 0.03 # the height of sample is 3cm
h_ana_center = 0.087 # unit m, distance from the center of analyzer sphere to the center of sample(i.e. origin)
h_array = np.linspace(h_ana_center-h_sample/2,h_ana_center+h_sample/2,31) # unit m, the distance from scattering point to spheric center of analyzer
R = 2.5 # unit m, radius of analyzer sphere
r_0 = 0.27 #unit m, radius of detector cylinder
h_d = 0.15 # unit m, the vertical offset of the detector center; elevation angle:15 cm -- 15.5 deg; 17.33 cm -- 20 deg;

# array assignment
gamma_deg = np.linspace(0,20,201)
gamma_rad = np.radians(gamma_deg)

for k in range(len(lambda_i_array)):
    lambda_i = lambda_i_array[k]
    name_Q = ['h_0']
    unit_Q = ['m']
    comment_Q = [f'{lambda_i*1e10}AA']
    data_Q = [h_0_inters]

    name_Q_reso = []
    unit_Q_reso = []
    comment_Q_reso = []
    data_Q_reso = []

    for j in range(len(theta_Dp_deg)):
        phi_f_deg = theta_Dp_deg[j]
        phi_f_rad = np.radians(phi_f_deg)   
        Q_centers = []
        Q_resos = []
        for i0 in range(len(h_0s)-1): # 对于每一个h_0s[i0]，我们都可以得到一个Q_center和Q_reso
            h_0_inter = (h_0s[i0+1]+h_0s[i0])/2
            Qs = []
            for i in range(len(h_array)):
                h = h_array[i]
                alpha_rad = sf.arcsin(h*np.cos(gamma_rad)/R)
                alpha_deg = np.degrees(alpha_rad)
                theta_B_deg = 90-alpha_deg
                theta_B_rad = np.radians(theta_B_deg)
                phi_deg = gamma_deg - alpha_deg + 90
                phi_rad = np.radians(phi_deg)
                beta_deg = 180 - phi_deg + alpha_deg
                beta_rad = np.radians(beta_deg)

                y_ana = h - R*np.cos(phi_rad)
                z_ana = R*np.sin(phi_rad)
                h_0 = y_ana - (z_ana - r_0)*sf.cot(beta_rad)
                y_off = h_0 - h_d
                rho = np.sqrt(r_0**2+(h_0-h)**2)

                psi_rad = sf.arctan((h_0-h)/r_0) 
                psi_deg = np.degrees(psi_rad)

                lambda_f = lambda_off*np.sin(theta_B_rad)
                theta_f_deg = 90 - gamma_deg
                theta_f_rad = np.radians(theta_f_deg)
                cos_2theta = np.sin(theta_i_rad)*np.sin(theta_f_rad)*np.cos(phi_i_rad-phi_f_rad)+np.cos(theta_i_rad)*np.cos(theta_f_rad)           
                twotheta_rad = sf.arccos(cos_2theta)
                twotheta_deg = np.degrees(twotheta_rad)
                Q = np.sqrt(4*np.pi**2*((1/(lambda_i**2))+(1/(lambda_f**2))-2*cos_2theta/(lambda_i*lambda_f)))

                for i1 in range(len(h_0)):
                    if h_0[i1]>=h_0s[i0] and h_0[i1]<=h_0s[i0+1]:
                        Qs.append(Q[i1])
            Q_min = min(Qs)
            Q_max = max(Qs)
            Q_center = (Q_max + Q_min)/2
            Q_reso = Q_max-Q_min
            Q_centers.append(Q_center)
            Q_resos.append(Q_reso)
        name_Q = name_Q+['Q',]
        unit_Q = unit_Q+['1/A',]
        comment_Q = comment_Q+[(f'{phi_f_deg}deg',)*1]
        data_Q = data_Q+[np.array(Q_centers)*1e-10]
        print (f'Loop-{k}{j} finished')
    # momentum output
    name_Q = _flatten(name_Q)
    unit_Q = _flatten(unit_Q)
    comment_Q = _flatten(comment_Q)
    data_Q = np.transpose(data_Q)
    filename = str(k+1).rjust(3,'0')+'_momentum transfer_'+f'lambda_i = {lambda_i*1e10} AA.dat'
    sf.dataoutput(name_Q,unit_Q,comment_Q,data_Q,datapath,filename)
print ('Run Successfully!')