From 58c02bd3b24cfbb399f1c68e6965f69b56e7e498 Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Thu, 18 Aug 2011 19:38:42 -0400 Subject: Replace Rayleigh scattering implementation with that from SNOMAN. The angular distribution is slightly different, and now fits with the distribution given in the GEANT4 physics reference manual. Unit test is now included to verify the correctness of the scattering. --- tests/rayleigh.py | 56 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 tests/rayleigh.py (limited to 'tests') diff --git a/tests/rayleigh.py b/tests/rayleigh.py new file mode 100644 index 0000000..4394ada --- /dev/null +++ b/tests/rayleigh.py @@ -0,0 +1,56 @@ +import unittest +import numpy as np + +from chroma.geometry import Solid, Geometry +from chroma.make import box +from chroma.sim import Simulation +from chroma.optics import water_wcsim +from chroma.event import Photons +import histogram +from histogram.root import rootify +import ROOT +ROOT.gROOT.SetBatch(1) + +class TestRayleigh(unittest.TestCase): + def setUp(self): + self.cube = Geometry() + self.cube.add_solid(Solid(box(100,100,100), water_wcsim, water_wcsim)) + self.cube.pmtids = [0] + self.sim = Simulation(self.cube, water_wcsim, bvh_bits=4, geant4_processes=0, + use_cache=False) + nphotons = 100000 + positions = np.tile([0,0,0], (nphotons,1)).astype(np.float32) + directions = np.tile([0,0,1], (nphotons,1)).astype(np.float32) + polarizations = np.zeros_like(positions) + phi = np.random.uniform(0, 2*np.pi, nphotons).astype(np.float32) + polarizations[:,0] = np.cos(phi) + polarizations[:,1] = np.sin(phi) + times = np.zeros(nphotons, dtype=np.float32) + wavelengths = np.empty(nphotons, np.float32) + wavelengths.fill(400.0) + + self.photons = Photons(positions=positions, directions=directions, polarizations=polarizations, + times=times, wavelengths=wavelengths) + + def testAngularDistributionPolarized(self): + # Fully polarized photons + self.photons.polarizations[:] = [1.0, 0.0, 0.0] + + photon_stop = self.sim.propagate_photons(self.photons, max_steps=1) + aborted = (photon_stop.histories & (1 << 31)) > 0 + self.assertFalse(aborted.any()) + + # Compute the dot product between initial and final directions + rayleigh_scatters = (photon_stop.histories & (1 << 4)) > 0 + cos_scatter = (self.photons.directions[rayleigh_scatters] * photon_stop.directions[rayleigh_scatters]).sum(axis=1) + theta_scatter = np.arccos(cos_scatter) + h = histogram.Histogram(bins=100, range=(0, np.pi)) + h.fill(theta_scatter) + h = rootify(h) + + # The functional form for polarized light should be (1 + \cos^2 \theta)\sin \theta + # according to GEANT4 physics reference manual + f = ROOT.TF1("pol_func", "[0]*(1+cos(x)**2)*sin(x)", 0, np.pi) + h.Fit(f) + self.assertGreater(f.GetProb(), 1e-3) + -- cgit From 6adecdaca7353b28e7a3d77ed563b5f98b86c137 Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Thu, 18 Aug 2011 19:39:43 -0400 Subject: Unit test to verify that photons at normal incidence do not abort. --- tests/propagation.py | 53 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100644 tests/propagation.py (limited to 'tests') diff --git a/tests/propagation.py b/tests/propagation.py new file mode 100644 index 0000000..331242b --- /dev/null +++ b/tests/propagation.py @@ -0,0 +1,53 @@ +import unittest +import numpy as np + +from chroma.geometry import Solid, Geometry +from chroma.make import box +from chroma.sim import Simulation +from chroma.optics import vacuum +from chroma.event import Photons + +class TestPropagation(unittest.TestCase): + def testAbort(self): + '''Photons that hit a triangle at normal incidence should not abort. + + Photons that hit a triangle at exactly normal incidence can sometimes produce a dot product + that is outside the range allowed by acos(). Trigger these with axis aligned photons in a box + ''' + + # Setup geometry + cube = Geometry() + cube.add_solid(Solid(box(100,100,100), vacuum, vacuum)) + cube.pmtids = [0] + sim = Simulation(cube, vacuum, bvh_bits=4, geant4_processes=0, + use_cache=False) + + # Create initial photons + nphotons = 10000 + positions = np.tile([0,0,0], (nphotons,1)).astype(np.float32) + directions = np.tile([0,0,1], (nphotons,1)).astype(np.float32) + polarizations = np.zeros_like(positions) + phi = np.random.uniform(0, 2*np.pi, nphotons).astype(np.float32) + polarizations[:,0] = np.cos(phi) + polarizations[:,1] = np.sin(phi) + times = np.zeros(nphotons, dtype=np.float32) + wavelengths = np.empty(nphotons, np.float32) + wavelengths.fill(400.0) + + photons = Photons(positions=positions, directions=directions, polarizations=polarizations, + times=times, wavelengths=wavelengths) + + # First make one step to check for strangeness + photon_stop = sim.propagate_photons(photons, max_steps=1) + self.assertFalse(np.isnan(photon_stop.positions).any()) + self.assertFalse(np.isnan(photon_stop.directions).any()) + self.assertFalse(np.isnan(photon_stop.polarizations).any()) + self.assertFalse(np.isnan(photon_stop.times).any()) + self.assertFalse(np.isnan(photon_stop.wavelengths).any()) + + # Now let it run the usual ten steps + photon_stop = sim.propagate_photons(photons, max_steps=10) + aborted = (photon_stop.histories & (1 << 31)) > 0 + print 'aborted photons: %1.1f' % (float(np.count_nonzero(aborted)) / nphotons) + self.assertFalse(aborted.any()) + -- cgit