import matplotlib.pyplot as plt
from numpy import sqrt, linspace, sum, random, array, concatenate, zeros_like, ones, ones_like, convolve
A = 0.275
E_0 = 25
c = 299792458 #m/s
m_e = 9.10938356e-31 #kg
E_e = linspace(0.1, E_0, 1000) 

def beta(E_e):
    return  A * E_e * sqrt(E_e**2 - (m_e * c**2)**2) * (E_0 - E_e)**2

def beta2(E_e, exponent=2.8):
    return  beta(E_e)**exponent

def plot_beta():
    fig, ax = plt.subplots()
    ax.plot(E_e, beta(E_e)/sum(beta(E_e)), c='k', lw=2, label=r'$\beta$ analytic')
    ax.set_xlabel(r'$E_e$', fontsize=14)
    ax.set_ylabel(r'$N(E_e)$', fontsize=14)
    ax.legend()
    plt.savefig('beta_decay.png', dpi=300, bbox_inches='tight')
    plt.show()
    
'''
Wir plotten ein analytisches Betaspektrum aenlich zu dem aus der Vorlesung 
auf Seite 49.
'''

plot_beta()

'''    
Beim 38Ca zerfall kommt es unter der Annahme des neutrinolosen 
doppelten Betazerfalls bei der registrierung beider Betaelektronen in einem 
Messintervall zu einem Deltapeak bei Q=(E_e1+E_e2)-2m_e im Energiespektrum.
Um dieses Energiespektrum experimentell nachzuweisen verwendet man Detektoren
und diese sind mit Messungenauigkeiten behaftet. 
Wir nehmen einen idealisierten Teilchendetektor mit sigma=5% Messungenauigkeit 
an um den Effekt zu visualisieren. Im wesentlichen stellt die Messung eine 
Faltung zwischen Deltafunktion und einer Gausseverteilung mit 
sigma=Energie*0.05 dar. 
'''

def measurement_uncertainty(energies, percent=5):
    return random.normal(energies, energies*percent/100)

'''
Um zunaechst ein Energiepektrum zweier gleichzeitig gemessenen Betaelektronen
darzustellen, die mit Dirac-Neutrinos $2/nu/beta/beta$ emittiert wurden und 
ohne eine analytische Representation der Gesamtenergie vorauszusetzen 
nutzen wir den Pseudozufallszahlengenerator von Python. So generieren wir eine 
Zufallsverteilung von einzelnen Betaelektronen.
Dies erleichtert es uns die Energieverteilung von gleichzeitig gemessenen 
Betaelektronen zu approximieren indem wir die zufellig generierten
Elektronenenergien addieren.
'''

def random_electron_energy(N):
    spectrum = beta(E_e)/sum(beta(E_e))
    electron_energies = []
    for i in list(range(len(spectrum)-1)):
        electron_energies.append(E_e[i] + random.rand(int(N * spectrum[i]))*(E_e[i+1]-E_e[i]))
    electron_energies = concatenate(array(electron_energies))
    random.shuffle(electron_energies)
    return electron_energies

electron_energies_1 = random_electron_energy(1e6)
electron_energies_2 = random_electron_energy(1e6)

E_0 = 50
E_e = linspace(0.1, E_0, 1000)

plt.figure()

plt.hist(electron_energies_1, bins=E_e, histtype='step', density=True, label=r'$\beta$ random')

plt.hist(electron_energies_2, bins=E_e, histtype='step', density=True, label=r'$\beta$ random')
double_beta_energies = electron_energies_1 + electron_energies_2

plt.hist(double_beta_energies,  bins=linspace(0.1, E_0, 1000), histtype='step', 
         density=True, label=r'$2\nu\beta\beta$')

E_0 = 25
E_e = linspace(0.1, E_0, 1000)

plt.plot(E_e, beta(E_e)/sum(beta(E_e))*40, c='k', lw=1, label=r'$\beta$ analytic')
plt.legend()
plt.savefig('double_beta_decay.png')
plt.show()

'''
Nun koennen wir die Funtkion fuer einzelne Elektronenzerfaelle im Exponenten
und E_0 anpassen, sodass sie mit dem Summenspektrum uebereinstimmt. Normalerweise wuerde 
man das nun mit einer Scipy-Funktion fitten und dabei die Residuen minimieren. Dies
wuerde den Umfang dieser Uebeung allerdings uebersteigen und wir geben uns mit 
einer Approximation der Approximation zufrieden.
'''

E_0 = 50
E_e = linspace(0.1, E_0, 1000)

plt.figure()
double_beta_energies = electron_energies_1 + electron_energies_2

plt.hist(double_beta_energies,  bins=linspace(0.1, E_0, 1000), histtype='step', 
         density=True, label=r'$2\nu\beta\beta$ random')

plt.plot(E_e, beta2(E_e)/sum(beta2(E_e))*20, c='r', lw=2, label=r'$2\nu\beta\beta$ analytic $N(E_e)^{2.8}$')

plt.legend()
plt.savefig('double_beta_decay_analytic.png')
plt.show()

'''
Wir sparen uns nun den Residuenplot, da wir wissen, dass es nicht praezise ist.
Dazu plotten wir den Deltapeak und verrauschen die Daten, so wie sie der Detektor
sehen wuerde.
'''

plt.figure()

double_beta_no_neutrino = zeros_like(E_e)
double_beta_no_neutrino[-1] = 0.1
electron_energies_no_neutrino = ones(100000)*E_0
weights_delta = ones_like(electron_energies_no_neutrino) * 0.000000001
weights = ones_like(electron_energies_no_neutrino) * 0.0000001

plt.hist(electron_energies_no_neutrino,  bins=linspace(0.1, E_0, 1000), histtype='step',weights=weights_delta, label=r'$0\nu\beta\beta$ random theoretical')

electron_energies_no_neutrino = measurement_uncertainty(electron_energies_no_neutrino)

plt.hist(electron_energies_no_neutrino,  bins=linspace(0.1, 2*E_0, 2000), histtype='step', weights=weights, label=r'$0\nu\beta\beta$ random measured')

double_beta_energies = electron_energies_1 + electron_energies_2

plt.hist(double_beta_energies,  bins=linspace(0.1, 2*E_0, 2000), histtype='step', 
         density=True, label=r'$2\nu\beta\beta$ random theoretical')

plt.hist(measurement_uncertainty(double_beta_energies),  bins=linspace(0.1, 2*E_0, 2000), histtype='step', 
         density=True, label=r'$2\nu\beta\beta$ random measured')

plt.plot(E_e, beta2(E_e)/sum(beta2(E_e))*20, c='r', lw=2, label=r'$2\nu\beta\beta$ analytic $N(E_e)^{2.8}$')

plt.xlim(40,55)
plt.ylim(0,0.001)
plt.legend()
plt.savefig('double_beta_decay_and_no_neutrino.png')
plt.show()