Demo entry 6718408

Back-Scattering

   

Submitted by anonymous on Mar 05, 2018 at 01:15
Language: Python 3. Code size: 2.0 kB.

import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

back_tally = 0
e_sum = 0
back_tally_square=0
e_sum_square=0

sims = 10000000
energies = np.zeros(sims)
energies = energies - 50
angles = np.zeros(sims)

for i in range(0,sims):
	random.seed(i)
	done = 0
	while done == 0:
		mu_0 = 2*random.uniform(0,1)-1

		rand_2 = random.uniform(0,1)

		e_prime = 1/(1+((1/0.511)*(1-mu_0)))
		e_ratio = e_prime/1

		ratio_KN_check = ((e_ratio**2)*((e_ratio**-1)+e_ratio-1+(mu_0**2)))/2

		if rand_2 <= ratio_KN_check:
			done = 1
			angles[i] = np.arccos(mu_0)
			if mu_0 >= -1 and mu_0 <= 0:
				back_tally += 1
				e_sum += e_prime
				back_tally_square += 1**2
				e_sum_square += e_prime**2
				energies[i] = e_prime
			

prob_back_scattered = back_tally/sims
ave_back_scattered_energy = e_sum/back_tally
RSD_back_scattered = (((sims/(sims-1))*(back_tally_square/(back_tally)**2))-((sims-1)**-1))**0.5
RSD_ave_back_scattered_energy = (((back_tally/(back_tally-1))*(e_sum_square/(e_sum)**2))-((back_tally-1)**-1))**0.5
print(prob_back_scattered)
print(ave_back_scattered_energy)
print(RSD_back_scattered)
print(RSD_ave_back_scattered_energy)

energies = energies[energies!=-50]

y, x = np.histogram(energies,bins=100,normed=True)
norm = np.sum(y)
y = y/norm
# integral=integrate.simps(y,x=x[:-1])
# print(integral)
# fig, ax = plt.subplots()
# ax.plot(x[:-1], y)
# plt.show()

def func(x, a, b, c):
	return a*np.exp(-b*x) + c

popt,pcov=curve_fit(func, x[:-1],  y)
print(popt[0],popt[1],popt[2])

plt.plot(x[:-1],y)
plt.plot(x[:-1],func(x[:-1],popt[0],popt[1],popt[2]))

equation=r'$y={%f}e^{-{%f}x}+{%f}$'%(popt[0],popt[1],popt[2])

plt.text(0.28, 1.5*np.average(y),equation, horizontalalignment='center',
     verticalalignment='center')
plt.title('Back-Scattering Energy Probability Distribution')
plt.xlabel('Energy (MeV)')
plt.ylabel('Relative Probability')
plt.show()

This snippet took 0.00 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).