diff options
Diffstat (limited to 'tests/rayleigh.py')
-rw-r--r-- | tests/rayleigh.py | 56 |
1 files changed, 0 insertions, 56 deletions
diff --git a/tests/rayleigh.py b/tests/rayleigh.py deleted file mode 100644 index 4394ada..0000000 --- a/tests/rayleigh.py +++ /dev/null @@ -1,56 +0,0 @@ -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) - |