python - python中的康普顿散射

Question

我和我的同事正在尝试使用 python 解决这个问题：

为了做到这一点，我们实现了这段代码：

#SAMPLING OF AN INTERACTION 

import numpy as np               #library for numerical calculations
import matplotlib.pyplot as plt  #library for plotting purposes
from scipy import constants      #needed for Physical and mathematical constants

#*******************************************************************************

c = 1                              #speed of light in natural units

me = constants.m_e*5.60958865*1e29  #electron mass in MeV/c^2

re = 2.8179403262e-13               #classical electron radius in cm   

E_init = 1                         #energy of incident photon in MeV

Na = constants.N_A                 #Avogadro number in 1/mol

d = 8.99*1e-5                      #hydrogen density in g/cm^3

iteration = 10**6                  #number of iteration  

A = 1                              #atomic weight of hydorgen in g/mol

N = (Na*d)/A                       #number of nuclei in 1/cm^3

t1 = np.random.rand(iteration)     #array of the given shape k and populate it 
                                   #with random samples from a uniform 
                                   #distribution over [0, 1).

t2=np.random.rand(iteration)       #array of the given shape d and populate it with 
                                   #random samples from a uniform distribution 
                                   #over [0, 1)

kappa = E_init/(me*c**2)           #generic position

theta=constants.pi*t2              #theta array over [0,pi]

E_f = E_init / (1+kappa*(1-np.cos(theta)))  #relation beetween theta and E_f

#*******************************************************************************

sigma_KN = (3/4)*(((1+kappa)/(kappa**3))*(((2*kappa*(1+kappa))/(1+2*kappa))-np.log(1+2*kappa)))+((1)/(2*kappa))*np.log(1+2*kappa)-((1+3*kappa)/(1+2*kappa)**2)  #total Klein-Nishima cross section in cm^2

mu = N*sigma_KN                    #interection lenght in cm^-1

s = - (1/mu)*np.log(t1)            #step lenght from interaction lenght

print("\nThe step lenght at", E_init, "MeV is:", s[0], "cm.\n" )
print("Meanwhile the step lenght distribution p(s) is: \n")

#*******************************************************************************

plt.figure()                       
plt.hist(s,bins=int(np.sqrt(iteration)), range=(0, 0.4*1E-18), density=True)
plt.grid()
plt.xlabel('s',fontweight='bold')
plt.ylabel('p(s)',fontweight='bold')
plt.title("Step lenght distribution")
plt.show()

print("\n")

#*******************************************************************************

E=np.linspace(E_init/(1+2*kappa),E_init,10000)

x=E_init/E

sigma_KN_diff=(constants.pi*(re**2)/(E_init*kappa))*((1/kappa)*(x-1)-1/x-x)

sigma_KN_total=(2*constants.pi*(re**2))*((((1+kappa)/(kappa**3))*(2*kappa*(1+kappa)/(1+2*kappa)-np.log(1+2*kappa)))+(1/(2*kappa))*np.log(1+2*kappa)-(1+3*kappa)/((1+2*kappa)**2))

sigma_KN_pdf=abs(sigma_KN_diff)/sigma_KN_total

print("The Klein-Nishima PDF is the following: \n")

#*******************************************************************************

plt.figure()
plt.plot(E,sigma_KN_pdf,'r--')
plt.xlabel('Final energy [MeV]')
plt.ylabel('Klein-Nishima PDF [$MeV^{-1}$]')
plt.grid()
plt.show()

#*******************************************************************************

E0=E[0]
E1=E[len(E)-1]
s0=sigma_KN_pdf[0]
s1=sigma_KN_pdf[len(sigma_KN_pdf)-1]

A=(1/2)*(s1+s0)*(E1-E0)
m=(s1-s0)/(E1-E0)
c=s0-E0*m
env_pdf=(1/A)*(m*E+c)
env_pdf_sc=1.24*env_pdf

print("\nThe envelop PDF of Klein-Nishima PDF is: \n")

#*******************************************************************************

plt.figure()
plt.plot(E,sigma_KN_pdf,'r--',label='Klein-Nishina PDF')
plt.plot(E,env_pdf_sc,'b--',label='Envelope PDF')
#plt.plot(E,env_pdf_sc-sigma_KN_pdf,'g--',label='Difference')
plt.xlabel('Final energy (MeV)')
plt.ylabel('Probability distribution ($MeV^{-1}$)')
plt.legend(loc='best')
plt.grid()
plt.show()

#*******************************************************************************
#Apply acceptance-rejection method

def sigma_fn(x):

    y=E_init/x

    return abs((constants.pi*(re**2)/(E_init*kappa))*((1/kappa)*(y-1)-1/y-y))/sigma_KN_total

def env_fn(x):

    return (1/A)*(m*x+c)

n=10**6

t1=np.random.rand(n)
t2=np.random.rand(n)

x=-np.sqrt(2*A*t1/m+((E0+c/m)**2))-(c/m)

x_sel=[]

for i in range(len(x)):
    if (t2[i]<=(sigma_fn(x[i])/(1.24*env_fn(x[i])))):
        x_sel+=[x[i]]

print("\nThe result of rejection Method is the following: \n")

#*******************************************************************************

plt.figure()
plt.hist(x_sel,bins=int(np.sqrt(n)),density=True)
plt.grid()
plt.xlabel('Final energy [MeV]')
plt.ylabel('Probability distribution [$MeV^{-1}$]')
plt.show()

#*******************************************************************************
#Convert final energy to scattering angle and plot

x_conv=np.divide(E_init,x_sel,dtype=float)
theta=np.arccos(1+(1/kappa)*(1-x_conv))

print("\nThe PDF of scattering angle of photon is: \n")

#*******************************************************************************

plt.figure()
plt.hist(theta,bins=int(np.sqrt(n)),density=True)
plt.grid()
plt.xlabel('Scattering angle [$\u03B8$]')
plt.ylabel('Probability distribution [$str^{-1}$]')
plt.show()

#*******************************************************************************
#Calculate electron properties

x1=np.array(x_sel,dtype=float)
theta1=np.array(theta,dtype=float)


E_ef = E_init-x1+me

phi = np.arctan(x1*np.sin(theta1)/(E_init-x1*np.cos(theta1)))                #final angle of electron

#E_ef = np.sqrt((E_init-x_sel[0]+me*c**2)**2-(me**2)*c**4)          #final energy of electron

#phi = np.arctan(1+kappa)*np.tan(theta1[0]/2)                #final angle of electron


print("\n The PDF of scattered electron is: \n")

#*******************************************************************************

plt.figure()
plt.hist(E_ef,bins=int(np.sqrt(n)),density=True)
plt.grid()
plt.xlabel('Final electron energy [MeV]')
plt.ylabel('Probability distribution [$MeV^{-1}$]')
plt.show()

print("\n The PDF of scattering angle of scattered electron is: \n")

plt.figure()
plt.hist(phi,bins=int(np.sqrt(n)),range=(0,3),density=True)
plt.grid()
plt.xlabel('Scattering angle [$\u03C6$]')
plt.ylabel('Probability distribution [$str^{-1}$]')
plt.show()

#*******************************************************************************

print("\nThe final state of scattered photon is: (",E_f[0]," MeV,", theta[0]," rad ) \n")



print("The final state of scattered electron is: (",E_ef[0]," MeV,", phi[0]," rad ) \n") # calculated from 4-momentum conservation (see https://en.wikipedia.org/wiki/Compton_scattering)

print("The energy conservation is verified?", E_init+me*c**2, "MeV is egual to" , x_sel[0]+E_ef[0], "MeV ? If", (E_init+me*c**2)/(x_sel[0]+E_ef[0]) ,"is 1 YES, otherwise NO. \n")

#*******************************************************************************

第一部分已完全解决（导出交互长度）。现在我们尝试使用差分横截面的拒绝方法对最终能量进行采样。

但是获得的结果（PDF VS E'）并不完全正确（它的行为众所周知）：

此外，能量守恒没有保留......

问题出在哪里？

标签： pythonphysics