From f6ec6d8fff619bdbc569e88cc48db20d74c2020e Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Fri, 19 Aug 2011 12:10:55 -0400 Subject: Rename test scripts so nose will find them. --- tests/propagation.py | 53 -------------------------------------------- tests/rayleigh.py | 56 ----------------------------------------------- tests/test_propagation.py | 53 ++++++++++++++++++++++++++++++++++++++++++++ tests/test_rayleigh.py | 56 +++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 109 insertions(+), 109 deletions(-) delete mode 100644 tests/propagation.py delete mode 100644 tests/rayleigh.py create mode 100644 tests/test_propagation.py create mode 100644 tests/test_rayleigh.py diff --git a/tests/propagation.py b/tests/propagation.py deleted file mode 100644 index 331242b..0000000 --- a/tests/propagation.py +++ /dev/null @@ -1,53 +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 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()) - 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) - diff --git a/tests/test_propagation.py b/tests/test_propagation.py new file mode 100644 index 0000000..331242b --- /dev/null +++ b/tests/test_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()) + diff --git a/tests/test_rayleigh.py b/tests/test_rayleigh.py new file mode 100644 index 0000000..4394ada --- /dev/null +++ b/tests/test_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 17f363a11109ed8ddeb983dcb6d47945ffb85f17 Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Fri, 19 Aug 2011 14:08:37 -0400 Subject: Liquid argon properties from MiniCLEAN --- detectors/miniclean.py | 4 ++-- optics.py | 13 +++++++++++++ 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/detectors/miniclean.py b/detectors/miniclean.py index 9223548..b0d4955 100644 --- a/detectors/miniclean.py +++ b/detectors/miniclean.py @@ -42,7 +42,7 @@ def build_miniclean(real_av=False): geo = Geometry() simple_iv = sphere(0.818) - geo.add_solid(Solid(simple_iv, water, vacuum, color=0x33FF0000)) + geo.add_solid(Solid(simple_iv, liquid_argon, vacuum, color=0x33FF0000)) polygons = read_polygons(os.path.join(dir, 'miniclean_polygons.txt')) @@ -69,7 +69,7 @@ def build_miniclean(real_av=False): geo.pmt_id_map = {} for i, cassette in enumerate(cassettes): polygon_mesh, polygon_colors = polygon_types[cassette['type']] - solid = Solid(polygon_mesh, water, water, surface=shiny_surface, + solid = Solid(polygon_mesh, liquid_argon, liquid_argon, surface=shiny_surface, color=polygon_colors) chroma_id = geo.add_solid(solid, cassette['rotation'], diff --git a/optics.py b/optics.py index 02b8e0d..7c89c08 100644 --- a/optics.py +++ b/optics.py @@ -242,3 +242,16 @@ water_wcsim.set('scattering_length', 1675.064, 1422.710, 1200.004, 1004.528, 833.9666, 686.1063])[::-1] / 100.0 * 0.625 # reversed, cm -> m, * magic tuning constant ) + +###### MiniCLEAN materials ###### + +liquid_argon = Material('liquid_argon') +liquid_argon.set('refractive_index', + wavelengths=[60.0,130.0,140.0,150.0,160.0,170.0,180.0,190.0,200.0,210.0,220.0,230.0,240.0,250.0,260.0,270.0,280.0,290.0,300.0,310.0,320.0,330.0,340.0,350.0,360.0,361.2,365,406.3,435.8,475.3,508.6,546.1,578,643.9,650,660,670,680,690,700,710,720,730,740,750,760,770,780,790,800], + value=[1.35247,1.35247,1.32727,1.30942,1.2962,1.28607,1.2781,1.27169,1.26645,1.2621,1.25845,1.25534,1.25267,1.25037,1.24835,1.24659,1.24503,1.24365,1.24242,1.24131,1.24032,1.23942,1.2386,1.23786,1.23719,1.237,1.2367,1.2347,1.2336,1.2324,1.2316,1.2308,1.2303,1.2295,1.22939,1.22929,1.22919,1.2291,1.22901,1.22893,1.22885,1.22877,1.2287,1.22863,1.22856,1.2285,1.22843,1.22837,1.22832,1.22826]) + +liquid_argon.set('absorption_length', 1e6) +liquid_argon.set('scattering_length', + wavelengths=[ 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0, 170.0, 180.0, 190.0, 200.0, 210.0, 220.0, 230.0, 240.0, 250.0, 260.0, 270.0, 280.0, 290.0, 300.0, 310.0, 320.0, 330.0, 340.0, 350.0, 360.0, 370.0, 380.0, 390.0, 400.0, 410.0, 420.0, 430.0, 440.0, 450.0, 800.0 ], + value=[ 43.5, 80.5, 137.3, 220.0, 335.3, 490.9, 695.2, 957.6, 1288.0, 1697.3, 2197.3, 2800.3, 3519.6, 4369.4, 5364.4, 6520.5, 7854.0, 9382.4, 11123.7, 13096.7, 15321.3, 17817.9, 20607.9, 23713.4, 27157.4, 30963.5, 35156.2, 39761.1, 44804.2, 50312.4, 56313.5, 62836.1, 69909.6, 77564.2, 85830.7, 94741.0, 104327.7, 114624.2, 125664.7, 137484.2, 1373291.0 ] + ) -- cgit