add simple bulk reemission

The ``Material`` struct now includes two new arrays: ``reemission_prob`` and ``reemission_cdf``. The former is sampled only when a photon is absorbed, and should be normalized accordingly. The latter defines the distribution from which the reemitted photon wavelength is drawn.

This process changes the photon wavelength in place, and is not capable of producing multiple secondaries. It also does not enforce energy conservation; the reemission spectrum is not itself wavelength-dependent.

1 files changed, 77 insertions, 0 deletions

diff --git a/test/test_reemission.py b/test/test_reemission.py
new file mode 100644
index 0000000..c6854be
--- /dev/null
+++ b/test/test_reemission.py
@@ -0,0 +1,77 @@
+import unittest
+import numpy as np
+import scipy.stats
+#import matplotlib.pyplot as plt
+
+from chroma.geometry import Solid, Geometry, Surface, Material
+from chroma.loader import create_geometry_from_obj
+from chroma.make import sphere
+from chroma.sim import Simulation
+from chroma.demo.optics import vacuum
+from chroma.event import Photons, SURFACE_DETECT
+from chroma.tools import enable_debug_on_crash
+
+class TestReemission(unittest.TestCase):
+    def testBulkReemission(self):
+        '''Test bulk reemission by starting a bunch of monoenergetic photons at
+        the center of a wavelength-shifting sphere, forcing reemission, and
+        checking that the final wavelength distribution matches the wls
+        spectrum.
+        '''
+        nphotons = 1000
+
+        # set up detector -- a sphere of 'scintillator' surrounded by a
+        # detecting sphere
+        scint = Material('scint')
+        scint.set('refractive_index', 1)
+        scint.set('absorption_length', 1)
+        scint.set('scattering_length', 1e7)
+        scint.set('reemission_prob', 1)
+
+        x = np.arange(0,1000,10)
+        norm = scipy.stats.norm(scale=50, loc=600)
+        pdf = nphotons * 10 * norm.pdf(x)
+        cdf = norm.cdf(x)
+        scint.reemission_cdf = np.array(zip(x, cdf))
+
+        detector = Surface('detector')
+        detector.set('detect', 1)
+
+        world = Geometry(vacuum)
+        world.add_solid(Solid(sphere(1000), vacuum, vacuum, surface=detector))
+        world.add_solid(Solid(sphere(100), scint, vacuum))
+        w = create_geometry_from_obj(world)
+
+        sim = Simulation(w, geant4_processes=0)
+
+        # initial photons -- isotropic 250 nm at the origin
+        pos = np.tile([0,0,0], (nphotons,1)).astype(np.float32)
+        dir = np.random.rand(nphotons, 3).astype(np.float32) * 2 - 1
+        dir /= np.sqrt(dir[:,0]**2 + dir[:,1]**2 + dir[:,2]**2)[:,np.newaxis]
+        pol = np.zeros_like(pos)
+        t = np.zeros(nphotons, dtype=np.float32)
+        wavelengths = np.ones(nphotons).astype(np.float32) * 250
+
+        photons = Photons(pos=pos, dir=dir, pol=pol, t=t, wavelengths=wavelengths)
+
+        # run simulation and extract final wavelengths
+        event = sim.simulate([photons], keep_photons_end=True).next()
+        mask = (event.photons_end.flags & SURFACE_DETECT) > 0
+        final_wavelengths = event.photons_end.wavelengths[mask]
+
+        # compare wavelength distribution to scintillator's reemission pdf
+        hist, edges = np.histogram(final_wavelengths, bins=x)
+        print 'detected', sum(hist), 'of', nphotons, 'photons'
+
+        chi2 = scipy.stats.chisquare(hist, pdf[:-1])[1]
+        print 'chi2 =', chi2
+
+        # show histogram comparison
+        #plt.figure(1)
+        #width = edges[1] - edges[0]
+        #plt.bar(left=edges, height=pdf, width=width, color='red')
+        #plt.bar(left=edges[:-1], height=hist, width=width)
+        #plt.show()
+
+        self.assertTrue(chi2 > 0.75)
+