Блог

Я и мой коллега пытаемся решить эту проблему с помощью 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’) не полностью корректен (его поведение хорошо известно):

Also, the energy conservation is not preserved…

Where is the problem?