diff options
| -rw-r--r-- | geometry.py | 32 | ||||
| -rw-r--r-- | pi0.py | 102 | ||||
| -rw-r--r-- | solids/pmts.py | 2 |
3 files changed, 114 insertions, 22 deletions
diff --git a/geometry.py b/geometry.py index ac71429..d80b730 100644 --- a/geometry.py +++ b/geometry.py @@ -49,8 +49,13 @@ class Mesh(object): self.vertices = np.array(unique_vertices) - def __getitem__(self, key): - return self.vertices[self.triangles[key]] + def assemble(self, key=slice(None), group=True): + if group: + vertex_indices = self.triangles[key] + else: + vertex_indices = self.triangles[key].flatten() + + return self.vertices[vertex_indices] def __len__(self): return len(self.triangles) @@ -167,15 +172,6 @@ def morton_order(mesh, bits): bits of the quantized center coordinates of each triangle. Each coordinate axis is quantized into 2**bits bins. """ - - #lower_bound = np.array([np.min(mesh[:,:,0]), - # np.min(mesh[:,:,1]), - # np.min(mesh[:,:,2])]) - - #upper_bound = np.array([np.max(mesh[:,:,0]), - # np.max(mesh[:,:,1]), - # np.max(mesh[:,:,2])]) - lower_bound, upper_bound = mesh.get_bounds() if bits <= 0 or bits > 10: @@ -186,8 +182,7 @@ def morton_order(mesh, bits): def quantize(x): return np.uint32((x-lower_bound)*max_value/(upper_bound-lower_bound)) - #mean_positions = quantize(np.mean(mesh, axis=1)) - mean_positions = quantize(np.mean(mesh[:], axis=1)) + mean_positions = quantize(np.mean(mesh.assemble(), axis=1)) return interleave(mean_positions, bits) @@ -288,14 +283,9 @@ class Geometry(object): i1 = np.searchsorted(zvalues_mesh, z) i2 = np.searchsorted(zvalues_mesh, z, side='right') - self.lower_bounds[i] = [np.min(self.mesh[i1:i2][:,:,0]), - np.min(self.mesh[i1:i2][:,:,1]), - np.min(self.mesh[i1:i2][:,:,2])] - - self.upper_bounds[i] = [np.max(self.mesh[i1:i2][:,:,0]), - np.max(self.mesh[i1:i2][:,:,1]), - np.max(self.mesh[i1:i2][:,:,2])] - + self.lower_bounds[i] = np.min(self.mesh.assemble(slice(i1,i2), group=False), axis=0) + self.upper_bounds[i] = np.max(self.mesh.assemble(slice(i1,i2), group=False), axis=0) + self.node_map[i] = i1 self.node_length[i] = i2-i1 self.layers[i] = layer @@ -0,0 +1,102 @@ +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 + +def rocket_to_lab(energy, momentum, v): + """ + Return the energy and momentum of a particle in the lab frame from its + energy and momentum in a rocket frame traveling at a velocity `v` with + 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. + """ + e0 = float(energy)#/c**2 + p0 = np.asarray(momentum, float)#/c + v = np.asarray(v, float)#/c + + try: + assert e0**2 - p0.dot(p0) >= -1.0e-70 + except AssertionError: + print e0**2 - p0.dot(p0) + raise + + g = 1.0/np.sqrt(1.0-v.dot(v)) + + x = np.dot(p0, v)/np.linalg.norm(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): + """ + Return the energy and directions for two photons produced from a pi0 + decay in the lab frame given that one of the photons decayed at polar + angles `theta` and `phi` in the pi0 rest frame. + + Args: + - energy: float + The total energy of the pi0 in MeV. + - direction: array + The direction of the pi0's velocity 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_v = pi0_p/pi0_e + + 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() diff --git a/solids/pmts.py b/solids/pmts.py index 377024d..22d9d67 100644 --- a/solids/pmts.py +++ b/solids/pmts.py @@ -51,7 +51,7 @@ def build_pmt(filename, glass_thickness, outer_material=water, theta=np.pi/8): outer_envelope = Solid(outer_envelope_mesh, glass, outer_material) - photocathode = np.mean(inner_envelope_mesh[:], axis=1)[:,1] > 0 + photocathode = np.mean(inner_envelope_mesh.assemble(), axis=1)[:,1] > 0 inner_envelope = Solid(inner_envelope_mesh, vacuum, glass, surface=np.where(photocathode, r7081hqe_photocathode, shiny_surface), color=np.where(photocathode, 0xff00, 0xff0000)) |