1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
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)
|
