summaryrefslogtreecommitdiff
path: root/tests/rayleigh.py
diff options
context:
space:
mode:
Diffstat (limited to 'tests/rayleigh.py')
-rw-r--r--tests/rayleigh.py56
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)
-