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.

6 files changed, 140 insertions, 43 deletions

diff --git a/chroma/cuda/geometry_types.h b/chroma/cuda/geometry_types.h
index 46226b2..6da1a47 100644
--- a/chroma/cuda/geometry_types.h
+++ b/chroma/cuda/geometry_types.h
@@ -6,17 +6,14 @@ struct Material
     float *refractive_index;
     float *absorption_length;
     float *scattering_length;
+    float *reemission_prob;
+    float *reemission_cdf;
     unsigned int n;
     float step;
     float wavelength0;
 };
 
-// surface models
-enum {
-    SURFACE_DEFAULT,  // specular + diffuse + absorption + detection
-    SURFACE_COMPLEX,  // use complex index of refraction
-    SURFACE_WLS       // wavelength-shifting reemission
-};
+enum { SURFACE_DEFAULT, SURFACE_COMPLEX, SURFACE_WLS };
 
 struct Surface
 {
diff --git a/chroma/cuda/hybrid_render.cu b/chroma/cuda/hybrid_render.cu
index 2c64d0e..29edefa 100644
--- a/chroma/cuda/hybrid_render.cu
+++ b/chroma/cuda/hybrid_render.cu
@@ -44,9 +44,6 @@ to_diffuse(Photon &p, State &s, Geometry *g, curandState &rng, int max_steps)
 	    if (p.history & REFLECT_DIFFUSE)
 		break;
 
-            if (p.history & SURFACE_REEMIT)
-                break;
-
 	    if (command == BREAK)
 		break;
 
diff --git a/chroma/cuda/photon.h b/chroma/cuda/photon.h
index 4ec1cea..cb1c092 100644
--- a/chroma/cuda/photon.h
+++ b/chroma/cuda/photon.h
@@ -36,6 +36,8 @@ struct State
     float refractive_index1, refractive_index2;
     float absorption_length;
     float scattering_length;
+    float reemission_prob;
+    Material *material1;
 
     int surface_index;
 
@@ -53,6 +55,7 @@ enum
     REFLECT_SPECULAR = 0x1 << 6,
     SURFACE_REEMIT   = 0x1 << 7,
     SURFACE_TRANSMIT = 0x1 << 8,
+    BULK_REEMIT      = 0x1 << 9,
     NAN_ABORT        = 0x1 << 31
 }; // processes
 
@@ -77,53 +80,56 @@ __device__ void
 fill_state(State &s, Photon &p, Geometry *g)
 {
     p.last_hit_triangle = intersect_mesh(p.position, p.direction, g,
-					 s.distance_to_boundary,
-					 p.last_hit_triangle);
+                                         s.distance_to_boundary,
+                                         p.last_hit_triangle);
 
     if (p.last_hit_triangle == -1) {
-	p.history |= NO_HIT;
-	return;
+        p.history |= NO_HIT;
+        return;
     }
-    
+
     Triangle t = get_triangle(g, p.last_hit_triangle);
-    
+
     unsigned int material_code = g->material_codes[p.last_hit_triangle];
-    
+
     int inner_material_index = convert(0xFF & (material_code >> 24));
     int outer_material_index = convert(0xFF & (material_code >> 16));
     s.surface_index = convert(0xFF & (material_code >> 8));
-    
+
     float3 v01 = t.v1 - t.v0;
     float3 v12 = t.v2 - t.v1;
-    
+
     s.surface_normal = normalize(cross(v01, v12));
-				 
+
     Material *material1, *material2;
     if (dot(s.surface_normal,-p.direction) > 0.0f) {
-	// outside to inside
-	material1 = g->materials[outer_material_index];
-	material2 = g->materials[inner_material_index];
+        // outside to inside
+        material1 = g->materials[outer_material_index];
+        material2 = g->materials[inner_material_index];
 
-	s.inside_to_outside = false;
+        s.inside_to_outside = false;
     }
     else {
-	// inside to outside
-	material1 = g->materials[inner_material_index];
-	material2 = g->materials[outer_material_index];
-	s.surface_normal = -s.surface_normal;
+        // inside to outside
+        material1 = g->materials[inner_material_index];
+        material2 = g->materials[outer_material_index];
+        s.surface_normal = -s.surface_normal;
 
-	s.inside_to_outside = true;
+        s.inside_to_outside = true;
     }
 
     s.refractive_index1 = interp_property(material1, p.wavelength,
-					  material1->refractive_index);
+                                          material1->refractive_index);
     s.refractive_index2 = interp_property(material2, p.wavelength,
-					  material2->refractive_index);
+                                          material2->refractive_index);
     s.absorption_length = interp_property(material1, p.wavelength,
-					  material1->absorption_length);
+                                          material1->absorption_length);
     s.scattering_length = interp_property(material1, p.wavelength,
-					  material1->scattering_length);
+                                          material1->scattering_length);
+    s.reemission_prob = interp_property(material1, p.wavelength,
+                                        material1->reemission_prob);
 
+    s.material1 = material1;
 } // fill_state
 
 __device__ float3
@@ -189,7 +195,7 @@ int propagate_to_boundary(Photon &p, State &s, curandState &rng,
     float absorption_distance = -s.absorption_length*logf(curand_uniform(&rng));
     float scattering_distance = -s.scattering_length*logf(curand_uniform(&rng));
 
-    if (use_weights && p.weight > WEIGHT_LOWER_THRESHOLD) // Prevent absorption
+    if (use_weights && p.weight > WEIGHT_LOWER_THRESHOLD && s.reemission_prob == 0) // Prevent absorption
 	absorption_distance = 1e30;
     else
 	use_weights = false;
@@ -224,15 +230,25 @@ int propagate_to_boundary(Photon &p, State &s, curandState &rng,
     }
 
     if (absorption_distance <= scattering_distance) {
-	if (absorption_distance <= s.distance_to_boundary) {
-	    p.time += absorption_distance/(SPEED_OF_LIGHT/s.refractive_index1);
-	    p.position += absorption_distance*p.direction;
-	    p.history |= BULK_ABSORB;
-
-	    p.last_hit_triangle = -1;
-
-	    return BREAK;
-	} // photon is absorbed in material1
+        if (absorption_distance <= s.distance_to_boundary) {
+            p.time += absorption_distance/(SPEED_OF_LIGHT/s.refractive_index1);
+            p.position += absorption_distance*p.direction;
+
+            float uniform_sample_reemit = curand_uniform(&rng);
+            if (uniform_sample_reemit < s.reemission_prob) {
+                p.wavelength = sample_cdf(&rng, s.material1->n, 
+                                          s.material1->wavelength0,
+                                          s.material1->step,
+                                          s.material1->reemission_cdf);
+                p.history |= BULK_REEMIT;
+                return CONTINUE;
+            } // photon is reemitted
+            else {
+                p.last_hit_triangle = -1;
+                p.history |= BULK_ABSORB;
+                return BREAK;
+            } // photon is absorbed in material1
+	}
     }
     else {
 	if (scattering_distance <= s.distance_to_boundary) {
@@ -257,7 +273,7 @@ int propagate_to_boundary(Photon &p, State &s, curandState &rng,
     // Scale weight by absorption probability along this distance
     if (use_weights)
 	p.weight *= expf(-s.distance_to_boundary/s.absorption_length);
-    
+
     p.position += s.distance_to_boundary*p.direction;
     p.time += s.distance_to_boundary/(SPEED_OF_LIGHT/s.refractive_index1);
 
diff --git a/chroma/geometry.py b/chroma/geometry.py
index 19294be..601dcf7 100644
--- a/chroma/geometry.py
+++ b/chroma/geometry.py
@@ -147,6 +147,8 @@ class Material(object):
         self.refractive_index = None
         self.absorption_length = None
         self.scattering_length = None
+        self.set('reemission_prob', 0)
+        self.set('reemission_cdf', 0)
         self.density = 0.0 # g/cm^3
         self.composition = {} # by mass
 
diff --git a/chroma/gpu/geometry.py b/chroma/gpu/geometry.py
index ff71b56..41d52a2 100644
--- a/chroma/gpu/geometry.py
+++ b/chroma/gpu/geometry.py
@@ -46,15 +46,23 @@ class GPUGeometry(object):
             absorption_length_gpu = ga.to_gpu(absorption_length)
             scattering_length = interp_material_property(wavelengths, material.scattering_length)
             scattering_length_gpu = ga.to_gpu(scattering_length)
+            reemission_prob = interp_material_property(wavelengths, material.reemission_prob)
+            reemission_prob_gpu = ga.to_gpu(reemission_prob)
+            reemission_cdf = interp_material_property(wavelengths, material.reemission_cdf)
+            reemission_cdf_gpu = ga.to_gpu(reemission_cdf)
 
             self.material_data.append(refractive_index_gpu)
             self.material_data.append(absorption_length_gpu)
             self.material_data.append(scattering_length_gpu)
+            self.material_data.append(reemission_prob)
+            self.material_data.append(reemission_cdf)
 
             material_gpu = \
                 make_gpu_struct(material_struct_size,
                                 [refractive_index_gpu, absorption_length_gpu,
                                  scattering_length_gpu,
+                                 reemission_prob_gpu,
+                                 reemission_cdf_gpu,
                                  np.uint32(len(wavelengths)),
                                  np.float32(wavelength_step),
                                  np.float32(wavelengths[0])])
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)
+