diff options
| -rw-r--r-- | chroma/cuda/geometry_types.h | 9 | ||||
| -rw-r--r-- | chroma/cuda/hybrid_render.cu | 3 | ||||
| -rw-r--r-- | chroma/cuda/photon.h | 84 | ||||
| -rw-r--r-- | chroma/geometry.py | 2 | ||||
| -rw-r--r-- | chroma/gpu/geometry.py | 8 | ||||
| -rw-r--r-- | test/test_reemission.py | 77 |
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) + |