diff options
| -rw-r--r-- | pi0.py | 89 |
1 files changed, 28 insertions, 61 deletions
@@ -1,10 +1,7 @@ import numpy as np -#c = 2.99792458e8 -#h = 6.6260689633e-34 -#joules_per_MeV = 1.602176487e-19/1e-6 -kg_per_MeV = 1.782661758e-36/1e-6 -pi0_mass = 134.9766*kg_per_MeV +_kg_per_MeV = 1.782661758e-36/1e-6 +_pi0_mass = 134.9766*_kg_per_MeV def rocket_to_lab(energy, momentum, v): """ @@ -13,23 +10,18 @@ def rocket_to_lab(energy, momentum, v): respect to the lab frame. Args: - - energy: float - The energy of the particle in the rocket frame in kg. - - momentum: array - The momentum of the particle in the rocket frame in kg. - - v: array - The velocity of the rocket frame as seen in the lab frame in units - of c. + - energy: float, kg + The energy of the particle in the rocket frame. + - momentum: array, kg + The momentum of the particle in the rocket frame. + - v: array, units of c + The velocity of the rocket frame as seen in the lab frame. """ - e0 = float(energy)#/c**2 - p0 = np.asarray(momentum, float)#/c - v = np.asarray(v, float)#/c + e0 = float(energy) + p0 = np.asarray(momentum, float) + v = np.asarray(v, float) - try: - assert e0**2 - p0.dot(p0) >= -1.0e-70 - except AssertionError: - print e0**2 - p0.dot(p0) - raise + assert e0**2 - p0.dot(p0) >= -1.0e-70 g = 1.0/np.sqrt(1.0-v.dot(v)) @@ -37,7 +29,6 @@ def rocket_to_lab(energy, momentum, v): p = p0 + ((g-1.0)*x + g*np.linalg.norm(v)*e0)*v/np.linalg.norm(v) e = np.sqrt(e0**2 - p0.dot(p0) + p.dot(p)) - #return e*c**2, p*c return e, p def pi0_decay(energy, direction, theta, phi): @@ -47,56 +38,32 @@ def pi0_decay(energy, direction, theta, phi): angles `theta` and `phi` in the pi0 rest frame. Args: - - energy: float - The total energy of the pi0 in MeV. + - energy: float, MeV + The total energy of the pi0. - direction: array The direction of the pi0's velocity in the lab frame. + + Returns: + - e1: float, MeV + The energy of the first photon in the lab frame. + - v1: array, + The direction of the first photon in the lab frame. + - e2: float, MeV + The energy of the second photon in the lab frame. + - v2: array, + The direction of the second photon in the lab frame. """ direction = np.asarray(direction)/np.linalg.norm(direction) - pi0_e = float(energy)*kg_per_MeV#/c**2 - pi0_p = np.sqrt(pi0_e**2-pi0_mass**2)*direction + pi0_e = float(energy)*_kg_per_MeV + pi0_p = np.sqrt(pi0_e**2-_pi0_mass**2)*direction pi0_v = pi0_p/pi0_e - photon_e0 = pi0_mass/2.0 + photon_e0 = _pi0_mass/2.0 photon_p0 = photon_e0*np.array([np.cos(theta)*np.sin(phi), np.sin(theta)*np.sin(phi), np.cos(phi)]) - #print photon_p0/np.linalg.norm(photon_p0) - e1, p1 = rocket_to_lab(photon_e0, photon_p0, pi0_v) v1 = p1/np.linalg.norm(p1) e2, p2 = rocket_to_lab(photon_e0, -photon_p0, pi0_v) v2 = p2/np.linalg.norm(p2) - return (e1/kg_per_MeV, v1), (e2/kg_per_MeV, v2) - -if __name__ == '__main__': - import sys - - npi0s = 10000 - - pi0_e = 300.0*kg_per_MeV - pi0_p = np.sqrt(pi0_e**2 - pi0_mass**2) - - print 'pi0 e: %f MeV' % (pi0_e/kg_per_MeV) - print 'pi0 p: %f MeV' % (pi0_p/kg_per_MeV) - - e, cos_theta = [], [] - for i, (theta, phi) in enumerate(zip(np.random.rand(npi0s), np.arccos(np.random.uniform(-1.0, 1.0, size=npi0s)))): - print '\r%i' % (i+1), - sys.stdout.flush() - - (e1, v1), (e2, v2) = pi0_decay(pi0_e/kg_per_MeV, [0,0,1], theta, phi) - e += [e1, e2] - cos_theta += [v1.dot(v2)] - - import matplotlib.pyplot as plt - - plt.figure() - plt.title('Energy Distribution of Photons from a %.0f MeV Pi0 Decay' % (pi0_e/kg_per_MeV)) - plt.hist(e, 100) - plt.xlabel('Energy (MeV)') - plt.figure() - plt.title('Opening Angle between Photons in Lab from a %.0f MeV Pi0 Decay' % (pi0_e/kg_per_MeV)) - plt.hist(cos_theta, 100) - plt.xlabel('Cos($\theta$)') - plt.show() + return (e1/_kg_per_MeV, v1), (e2/_kg_per_MeV, v2) |