From 01d5527db39520ca548518ed1194f8b863a4f077 Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Wed, 10 Aug 2011 11:25:39 -0400 Subject: Rename chroma.io to chroma.fileio to avoid collision with Python package named io --- .rootlogon.C | 2 +- fileio/__init__.py | 0 fileio/root.C | 95 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ fileio/root.py | 19 +++++++++++ io/__init__.py | 0 io/root.C | 95 ------------------------------------------------------ io/root.py | 19 ----------- sim.py | 2 +- 8 files changed, 116 insertions(+), 116 deletions(-) create mode 100644 fileio/__init__.py create mode 100644 fileio/root.C create mode 100644 fileio/root.py delete mode 100644 io/__init__.py delete mode 100644 io/root.C delete mode 100644 io/root.py diff --git a/.rootlogon.C b/.rootlogon.C index f1b4cf6..17915c5 100644 --- a/.rootlogon.C +++ b/.rootlogon.C @@ -1,3 +1,3 @@ { - gROOT->ProcessLine(".L io/root.C+g"); + gROOT->ProcessLine(".L fileio/root.C+g"); } diff --git a/fileio/__init__.py b/fileio/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/fileio/root.C b/fileio/root.C new file mode 100644 index 0000000..096811a --- /dev/null +++ b/fileio/root.C @@ -0,0 +1,95 @@ +#include +#include +#include +#include + +struct Photon { + double t; + TVector3 pos; + TVector3 dir; + TVector3 pol; + double wavelength; // nm + unsigned int history; + int last_hit_triangle; +}; + +struct MC { + std::string particle; + TVector3 gen_pos; + TVector3 gen_dir; + double gen_total_e; + + std::vector photon_start; + std::vector photon_stop; + +}; + +struct Channel { + Channel() : channel_id(-1), t(-9999.0), q(-9999.0) { }; + int channel_id; + double t; + double q; + unsigned int mc_history; +}; + +struct Event { + int event_id; + MC mc; + int nhit; + std::vector channel; + +}; + +void fill_photons(Event *ev, bool start, + unsigned int nphotons, float *pos, float *dir, + float *pol, float *wavelength, float *t0, + int *histories=0, int *last_hit_triangle=0) +{ + std::vector &photons = start ? ev->mc.photon_start : ev->mc.photon_stop; + photons.resize(nphotons); + + for (unsigned int i=0; i < nphotons; i++) { + Photon &photon = photons[i]; + photon.t = t0[i]; + photon.pos.SetXYZ(pos[3*i], pos[3*i + 1], pos[3*i + 2]); + photon.dir.SetXYZ(dir[3*i], dir[3*i + 1], dir[3*i + 2]); + photon.pol.SetXYZ(pol[3*i], pol[3*i + 1], pol[3*i + 2]); + photon.wavelength = wavelength[i]; + if (histories) + photon.history = histories[i]; + else + photon.history = 0; + + if (last_hit_triangle) + photon.last_hit_triangle = last_hit_triangle[i]; + else + photon.last_hit_triangle = -1; + } +} + + +void fill_hits(Event *ev, unsigned int nchannels, float *t, + float *q, unsigned int *history) +{ + ev->channel.resize(0); + ev->nhit = 0; + Channel ch; + for (unsigned int i=0; i < nchannels; i++) { + if (t[i] < 1e8) { + ev->nhit++; + ch.channel_id = i; + ch.t = t[i]; + ch.q = q[i]; + ch.mc_history = history[i]; + ev->channel.push_back(ch); + } + } +} + + +#ifdef __MAKECINT__ +#pragma link C++ class vector; +#pragma link C++ class vector; +#endif + + diff --git a/fileio/root.py b/fileio/root.py new file mode 100644 index 0000000..e43f5d4 --- /dev/null +++ b/fileio/root.py @@ -0,0 +1,19 @@ +import ROOT +import os.path + +ROOT.gROOT.ProcessLine('.L '+os.path.join(os.path.dirname(__file__), 'root.C+g')) + +from ROOT import Event + +fill_photons = ROOT.fill_photons +fill_hits = ROOT.fill_hits + +def make_tree(name, desc=''): + '''Create a ROOT tree for holding event information. + + Returns tuple of Event object for filling and TTree. + ''' + tree = ROOT.TTree(name, desc) + ev = ROOT.Event() + tree.Branch('ev', ev) + return ev, tree diff --git a/io/__init__.py b/io/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/io/root.C b/io/root.C deleted file mode 100644 index 096811a..0000000 --- a/io/root.C +++ /dev/null @@ -1,95 +0,0 @@ -#include -#include -#include -#include - -struct Photon { - double t; - TVector3 pos; - TVector3 dir; - TVector3 pol; - double wavelength; // nm - unsigned int history; - int last_hit_triangle; -}; - -struct MC { - std::string particle; - TVector3 gen_pos; - TVector3 gen_dir; - double gen_total_e; - - std::vector photon_start; - std::vector photon_stop; - -}; - -struct Channel { - Channel() : channel_id(-1), t(-9999.0), q(-9999.0) { }; - int channel_id; - double t; - double q; - unsigned int mc_history; -}; - -struct Event { - int event_id; - MC mc; - int nhit; - std::vector channel; - -}; - -void fill_photons(Event *ev, bool start, - unsigned int nphotons, float *pos, float *dir, - float *pol, float *wavelength, float *t0, - int *histories=0, int *last_hit_triangle=0) -{ - std::vector &photons = start ? ev->mc.photon_start : ev->mc.photon_stop; - photons.resize(nphotons); - - for (unsigned int i=0; i < nphotons; i++) { - Photon &photon = photons[i]; - photon.t = t0[i]; - photon.pos.SetXYZ(pos[3*i], pos[3*i + 1], pos[3*i + 2]); - photon.dir.SetXYZ(dir[3*i], dir[3*i + 1], dir[3*i + 2]); - photon.pol.SetXYZ(pol[3*i], pol[3*i + 1], pol[3*i + 2]); - photon.wavelength = wavelength[i]; - if (histories) - photon.history = histories[i]; - else - photon.history = 0; - - if (last_hit_triangle) - photon.last_hit_triangle = last_hit_triangle[i]; - else - photon.last_hit_triangle = -1; - } -} - - -void fill_hits(Event *ev, unsigned int nchannels, float *t, - float *q, unsigned int *history) -{ - ev->channel.resize(0); - ev->nhit = 0; - Channel ch; - for (unsigned int i=0; i < nchannels; i++) { - if (t[i] < 1e8) { - ev->nhit++; - ch.channel_id = i; - ch.t = t[i]; - ch.q = q[i]; - ch.mc_history = history[i]; - ev->channel.push_back(ch); - } - } -} - - -#ifdef __MAKECINT__ -#pragma link C++ class vector; -#pragma link C++ class vector; -#endif - - diff --git a/io/root.py b/io/root.py deleted file mode 100644 index e43f5d4..0000000 --- a/io/root.py +++ /dev/null @@ -1,19 +0,0 @@ -import ROOT -import os.path - -ROOT.gROOT.ProcessLine('.L '+os.path.join(os.path.dirname(__file__), 'root.C+g')) - -from ROOT import Event - -fill_photons = ROOT.fill_photons -fill_hits = ROOT.fill_hits - -def make_tree(name, desc=''): - '''Create a ROOT tree for holding event information. - - Returns tuple of Event object for filling and TTree. - ''' - tree = ROOT.TTree(name, desc) - ev = ROOT.Event() - tree.Branch('ev', ev) - return ev, tree diff --git a/sim.py b/sim.py index 34477f9..7644266 100755 --- a/sim.py +++ b/sim.py @@ -8,7 +8,7 @@ import detectors import optics import gpu import g4gen -from io import root +from fileio import root import numpy as np import math import ROOT -- cgit From ed5b635e31464585aebae8f462c1369224c2681d Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Wed, 10 Aug 2011 12:40:45 -0400 Subject: Using gzip compression level 1, the BVH at 10 bit is 7x smaller. Adds 10 seconds to cache load time if the BVH file is in the memory cache, otherwise this is faster if the BVH file must be read fresh from disk. --- geometry.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/geometry.py b/geometry.py index 579e01d..35da2ac 100644 --- a/geometry.py +++ b/geometry.py @@ -5,6 +5,7 @@ from itertoolset import * from tools import timeit from hashlib import md5 import cPickle as pickle +import gzip # all material/surface properties are interpolated at these # wavelengths when they are sent to the gpu @@ -285,7 +286,7 @@ class Geometry(object): cache_path = os.path.join(cache_dir, cache_file) try: - f = open(cache_path, 'rb') + f = gzip.GzipFile(cache_path, 'rb') except IOError: pass else: @@ -357,13 +358,15 @@ class Geometry(object): if unique_zvalues.size == 1: break + print >>sys.stderr, 'Writing BVH to cache directory...' + if not os.path.exists(cache_dir): os.makedirs(cache_dir) - f = open(cache_path, 'wb') - data = {} - for key in ['material1_index', 'material2_index', 'surface_index', 'colors', 'solid_id', 'lower_bounds', 'upper_bounds', 'node_map', 'node_map_end', 'layers', 'first_node']: - data[key] = getattr(self, key) - data['reorder'] = reorder - pickle.dump(data, f, -1) - f.close() + with gzip.GzipFile(cache_path, 'wb', compresslevel=1) as f: + data = {} + for key in ['material1_index', 'material2_index', 'surface_index', 'colors', 'solid_id', 'lower_bounds', 'upper_bounds', 'node_map', 'node_map_end', 'layers', 'first_node']: + data[key] = getattr(self, key) + data['reorder'] = reorder + pickle.dump(data, f, -1) + -- cgit From 0dcfffe29e8e8baac6f854cb6fde0e09d81ab62a Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Wed, 10 Aug 2011 19:12:45 -0400 Subject: Fix misuse of Material.set(), so now wavelengths and values of material properties are not transposed. --- optics.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/optics.py b/optics.py index c1fcf0a..d65c5f9 100644 --- a/optics.py +++ b/optics.py @@ -157,14 +157,14 @@ glass.set('scattering_length', 1e6) # From SNO+ database acrylic_sno = Material('acrylic_sno') -acrylic_sno.set('refractive_index', [ 200.0, 220.0, 240.0, 260.0, 280.0, 300.0, 320.0, 340.0, 360.0, 380.0, 400.0, 420.0, 440.0, 460.0, 480.0, 500.0, 520.0, 540.0, 560.0, 580.0, 600.0, 620.0, 640.0, 660.0, 680.0, 700.0, 720.0, 740.0, 760.0, 780.0, 800.0 ], [ 1.59816, 1.57399, 1.55692, 1.54436, 1.5348, 1.52732, 1.52133, 1.51644, 1.51238, 1.50897, 1.50607, 1.50357, 1.50139, 1.49949, 1.49781, 1.49632, 1.49498, 1.49378, 1.49269, 1.49171, 1.49081, 1.48999, 1.48924, 1.48854, 1.4879, 1.4873, 1.48675, 1.48624, 1.48576, 1.48531, 1.48488 ]) -acrylic_sno.set('absorption_length', [ 200.0, 220.0, 240.0, 260.0, 280.0, 300.0, 320.0, 340.0, 360.0, 380.0, 400.0, 420.0, 440.0, 460.0, 480.0, 500.0, 520.0, 540.0, 560.0, 580.0, 600.0, 620.0, 640.0, 660.0, 680.0, 700.0, 720.0, 740.0, 760.0, 780.0, 800.0 ], [ 0.0456452, 0.0512064, 0.0583106, 0.0677037, 0.080704, 0.0998834, 0.131021, 0.190364, 0.347969, 0.661751, 0.979776, 1.31804, 1.34242, 1.36771, 1.39397, 1.42126, 1.42019, 1.41911, 1.41804, 1.41697, 1.4159, 1.41483, 1.41376, 1.4127, 1.41163, 1.41057, 1.40951, 1.40845, 1.40739, 1.40634, 1.40528 ]) -acrylic_sno.set('scattering_length', [ 200.0, 220.0, 240.0, 260.0, 280.0, 300.0, 320.0, 340.0, 360.0, 380.0, 400.0, 420.0, 440.0, 460.0, 480.0, 500.0, 520.0, 540.0, 560.0, 580.0, 600.0, 620.0, 640.0, 660.0, 680.0, 700.0, 720.0, 740.0, 760.0, 780.0, 800.0 ], [ 4.60015, 6.73532, 9.53941, 13.1394, 17.6734, 23.2903, 30.1503, 38.4246, 48.2953, 59.9557, 73.61, 89.4735, 107.773, 128.745, 152.638, 179.713, 210.238, 244.497, 282.781, 325.395, 372.652, 424.879, 482.413, 545.601, 614.801, 690.385, 772.733, 862.236, 959.298, 1064.33, 1177.77 ]) +acrylic_sno.set('refractive_index', wavelengths=[ 200.0, 220.0, 240.0, 260.0, 280.0, 300.0, 320.0, 340.0, 360.0, 380.0, 400.0, 420.0, 440.0, 460.0, 480.0, 500.0, 520.0, 540.0, 560.0, 580.0, 600.0, 620.0, 640.0, 660.0, 680.0, 700.0, 720.0, 740.0, 760.0, 780.0, 800.0 ], value=[ 1.59816, 1.57399, 1.55692, 1.54436, 1.5348, 1.52732, 1.52133, 1.51644, 1.51238, 1.50897, 1.50607, 1.50357, 1.50139, 1.49949, 1.49781, 1.49632, 1.49498, 1.49378, 1.49269, 1.49171, 1.49081, 1.48999, 1.48924, 1.48854, 1.4879, 1.4873, 1.48675, 1.48624, 1.48576, 1.48531, 1.48488 ]) +acrylic_sno.set('absorption_length', wavelengths=[ 200.0, 220.0, 240.0, 260.0, 280.0, 300.0, 320.0, 340.0, 360.0, 380.0, 400.0, 420.0, 440.0, 460.0, 480.0, 500.0, 520.0, 540.0, 560.0, 580.0, 600.0, 620.0, 640.0, 660.0, 680.0, 700.0, 720.0, 740.0, 760.0, 780.0, 800.0 ], value=[ 0.0456452, 0.0512064, 0.0583106, 0.0677037, 0.080704, 0.0998834, 0.131021, 0.190364, 0.347969, 0.661751, 0.979776, 1.31804, 1.34242, 1.36771, 1.39397, 1.42126, 1.42019, 1.41911, 1.41804, 1.41697, 1.4159, 1.41483, 1.41376, 1.4127, 1.41163, 1.41057, 1.40951, 1.40845, 1.40739, 1.40634, 1.40528 ]) +acrylic_sno.set('scattering_length', wavelengths=[ 200.0, 220.0, 240.0, 260.0, 280.0, 300.0, 320.0, 340.0, 360.0, 380.0, 400.0, 420.0, 440.0, 460.0, 480.0, 500.0, 520.0, 540.0, 560.0, 580.0, 600.0, 620.0, 640.0, 660.0, 680.0, 700.0, 720.0, 740.0, 760.0, 780.0, 800.0 ], value=[ 4.60015, 6.73532, 9.53941, 13.1394, 17.6734, 23.2903, 30.1503, 38.4246, 48.2953, 59.9557, 73.61, 89.4735, 107.773, 128.745, 152.638, 179.713, 210.238, 244.497, 282.781, 325.395, 372.652, 424.879, 482.413, 545.601, 614.801, 690.385, 772.733, 862.236, 959.298, 1064.33, 1177.77 ]) # From SNO+ database labppo_scintillator = Material('labppo_scintillator') -labppo_scintillator.set('refractive_index', [ 200.0, 213.0, 215.0, 217.0, 219.0, 221.0, 223.0, 225.0, 227.0, 229.0, 231.0, 233.0, 235.0, 237.0, 239.0, 241.0, 243.0, 245.0, 247.0, 249.0, 251.0, 253.0, 255.0, 257.0, 259.0, 261.0, 263.0, 265.0, 267.0, 269.0, 271.0, 273.0, 275.0, 277.0, 279.0, 281.0, 283.0, 285.0, 287.0, 289.0, 291.0, 293.0, 295.0, 297.0, 299.0, 301.0, 303.0, 305.0, 307.0, 309.0, 311.0, 313.0, 315.0, 317.0, 319.0, 321.0, 323.0, 325.0, 327.0, 329.0, 331.0, 333.0, 335.0, 337.0, 339.0, 341.0, 343.0, 345.0, 347.0, 349.0, 351.0, 353.0, 355.0, 357.0, 359.0, 361.0, 363.0, 365.0, 367.0, 369.0, 371.0, 373.0, 375.0, 377.0, 379.0, 381.0, 383.0, 385.0, 387.0, 389.0, 391.0, 393.0, 395.0, 397.0, 399.0, 401.0, 403.0, 405.0, 407.0, 409.0, 411.0, 413.0, 415.0, 417.0, 419.0, 421.0, 423.0, 425.0, 427.0, 429.0, 431.0, 433.0, 435.0, 437.0, 439.0, 441.0, 443.0, 445.0, 447.0, 449.0, 451.0, 453.0, 455.0, 457.0, 459.0, 461.0, 463.0, 465.0, 467.0, 469.0, 471.0, 473.0, 475.0, 477.0, 479.0, 481.0, 483.0, 485.0, 487.0, 489.0, 491.0, 493.0, 495.0, 497.0, 499.0, 501.0, 503.0, 505.0, 507.0, 509.0, 511.0, 513.0, 515.0, 517.0, 519.0, 521.0, 523.0, 525.0, 527.0, 529.0, 531.0, 533.0, 535.0, 537.0, 539.0, 541.0, 543.0, 545.0, 547.0, 549.0, 551.0, 553.0, 555.0, 557.0, 559.0, 561.0, 563.0, 565.0, 567.0, 569.0, 571.0, 573.0, 575.0, 577.0, 579.0, 581.0, 583.0, 585.0, 587.0, 589.0, 591.0, 593.0, 595.0, 597.0, 599.0, 601.0, 603.0, 605.0, 607.0, 609.0, 611.0, 613.0, 615.0, 617.0, 619.0, 621.0, 623.0, 625.0, 627.0, 629.0, 631.0, 633.0, 635.0, 637.0, 639.0, 641.0, 643.0, 645.0, 647.0, 649.0, 651.0, 653.0, 655.0, 657.0, 659.0, 661.0, 663.0, 665.0, 667.0, 669.0, 671.0, 673.0, 675.0, 677.0, 679.0, 681.0, 683.0, 685.0, 687.0, 689.0, 691.0, 693.0, 695.0, 697.0, 699.0, 701.0, 703.0, 705.0, 707.0, 709.0, 711.0, 713.0, 715.0, 717.0, 719.0, 721.0, 723.0, 725.0, 727.0, 729.0, 731.0, 733.0, 735.0, 737.0, 739.0, 741.0, 743.0, 745.0, 747.0, 749.0, 751.0, 753.0, 755.0, 757.0, 759.0, 761.0, 763.0, 765.0, 767.0, 769.0, 771.0, 773.0, 775.0, 777.0, 779.0, 781.0, 783.0, 785.0, 787.0, 789.0, 791.0, 793.0, 795.0, 797.0, 799.0, 800.0 ], [ 1.75541, 1.75541, 1.75541, 1.7526, 1.74552, 1.73503, 1.72202, 1.70736, 1.69194, 1.67662, 1.66228, 1.65227, 1.64307, 1.63526, 1.62846, 1.62243, 1.61699, 1.61204, 1.60748, 1.60326, 1.59933, 1.59565, 1.59218, 1.58892, 1.58582, 1.58288, 1.58009, 1.57742, 1.57488, 1.57244, 1.57011, 1.56787, 1.56571, 1.56365, 1.56165, 1.55973, 1.55788, 1.55609, 1.55437, 1.5527, 1.55108, 1.54952, 1.548, 1.54654, 1.54511, 1.54373, 1.54239, 1.54109, 1.53982, 1.53859, 1.5374, 1.53623, 1.5351, 1.534, 1.53293, 1.53188, 1.53086, 1.52987, 1.5289, 1.52795, 1.52703, 1.52613, 1.52525, 1.52439, 1.52355, 1.52273, 1.52193, 1.52115, 1.52039, 1.51964, 1.51891, 1.51819, 1.51749, 1.5168, 1.51613, 1.51547, 1.51483, 1.5142, 1.51358, 1.51297, 1.51238, 1.5118, 1.51122, 1.51066, 1.51012, 1.50958, 1.50905, 1.50853, 1.50802, 1.50752, 1.50703, 1.50655, 1.50608, 1.50561, 1.50516, 1.50471, 1.50427, 1.50383, 1.50341, 1.50299, 1.50258, 1.50218, 1.50178, 1.50139, 1.501, 1.50063, 1.50025, 1.49989, 1.49953, 1.49918, 1.49883, 1.49848, 1.49815, 1.49781, 1.49749, 1.49717, 1.49685, 1.49654, 1.49623, 1.49593, 1.49563, 1.49533, 1.49504, 1.49476, 1.49448, 1.4942, 1.49392, 1.49366, 1.49339, 1.49313, 1.49287, 1.49261, 1.49236, 1.49212, 1.49187, 1.49163, 1.49139, 1.49116, 1.49093, 1.4907, 1.49047, 1.49025, 1.49003, 1.48982, 1.4896, 1.48939, 1.48918, 1.48898, 1.48877, 1.48857, 1.48838, 1.48818, 1.48799, 1.4878, 1.48761, 1.48742, 1.48724, 1.48706, 1.48688, 1.4867, 1.48653, 1.48635, 1.48618, 1.48601, 1.48585, 1.48568, 1.48552, 1.48536, 1.4852, 1.48504, 1.48488, 1.48473, 1.48458, 1.48443, 1.48428, 1.48413, 1.48398, 1.48384, 1.4837, 1.48356, 1.48342, 1.48328, 1.48314, 1.48301, 1.48287, 1.48274, 1.48261, 1.48248, 1.48235, 1.48223, 1.4821, 1.48198, 1.48185, 1.48173, 1.48161, 1.48149, 1.48137, 1.48126, 1.48114, 1.48103, 1.48091, 1.4808, 1.48069, 1.48058, 1.48047, 1.48036, 1.48025, 1.48015, 1.48004, 1.47994, 1.47984, 1.47973, 1.47963, 1.47953, 1.47943, 1.47933, 1.47924, 1.47914, 1.47904, 1.47895, 1.47885, 1.47876, 1.47867, 1.47858, 1.47849, 1.4784, 1.47831, 1.47822, 1.47813, 1.47804, 1.47796, 1.47787, 1.47779, 1.4777, 1.47762, 1.47754, 1.47745, 1.47737, 1.47729, 1.47721, 1.47713, 1.47705, 1.47698, 1.4769, 1.47682, 1.47675, 1.47667, 1.4766, 1.47652, 1.47645, 1.47637, 1.4763, 1.47623, 1.47616, 1.47609, 1.47602, 1.47595, 1.47588, 1.47581, 1.47574, 1.47567, 1.4756, 1.47554, 1.47547, 1.4754, 1.47534, 1.47527, 1.47521, 1.47515, 1.47508, 1.47502, 1.47496, 1.47489, 1.47483, 1.47477, 1.47471, 1.47465, 1.47459, 1.47453, 1.47447, 1.47441, 1.47435, 1.47429, 1.47424, 1.47418, 1.47412, 1.47407, 1.47401, 1.47395, 1.4739, 1.47384, 1.47379, 1.47373, 1.47368, 1.47363, 1.4736 ]) -labppo_scintillator.set('absorption_length', [ 200.0, 215.0, 225.0, 235.0, 245.0, 255.0, 265.0, 275.0, 285.0, 295.0, 305.0, 315.0, 325.0, 335.0, 345.0, 355.0, 365.0, 375.0, 385.0, 395.0, 405.0, 415.0, 425.0, 435.0, 445.0, 455.0, 465.0, 475.0, 485.0, 495.0, 505.0, 515.0, 525.0, 535.0, 545.0, 555.0, 565.0, 575.0, 585.0, 595.0, 800.0 ], [ 1.54944e-07, 2.99637e-07, 1.72048e-05, 3.41222e-05, 1.53692e-05, 7.92675e-06, 1.00564e-05, 0.000270688, 0.000949206, 0.00134987, 0.00322232, 0.011032, 0.0379538, 0.708042, 2.68045, 4.20713, 5.71163, 7.92686, 9.90035, 17.239, 28.2318, 41.8939, 67.3786, 120.517, 347.953, 412.636, 147.016, 80.7207, 54.7561, 45.1483, 39.4147, 39.1438, 42.199, 41.6337, 35.783, 47.614, 58.3934, 43.0845, 32.0874, 35.6067, 35.6067 ]) -labppo_scintillator.set('absorption_wls_length', [ 200.0, 215.0, 225.0, 235.0, 245.0, 255.0, 265.0, 275.0, 285.0, 295.0, 305.0, 315.0, 325.0, 335.0, 345.0, 355.0, 365.0, 375.0, 385.0, 395.0, 405.0, 415.0, 425.0, 435.0, 445.0, 455.0, 465.0, 475.0, 485.0, 495.0, 505.0, 515.0, 525.0, 535.0, 545.0, 555.0, 565.0, 575.0, 585.0, 595.0, 800.0 ], [ 0.00109808, 0.000693615, 0.000437277, 0.000276096, 0.000173933, 0.000103044, 6.03065e-05, 3.56199e-05, 2.25651e-05, 1.60844e-05, 1.62228e-05, 2.21832e-05, 3.03067e-05, 5.24247e-05, 0.000484492, 0.0127541, 0.294651, 4.15892, 15.6626, 25.4821, 36.1763, 50.723, 63.9962, 77.7386, 89.8611, 103.551, 115.859, 125.442, 146.664, 163.748, 194.986, 194.652, 202.262, 234.033, 262.861, 279.883, 323.991, 299.632, 287.367, 303.42, 303.42 ]) -labppo_scintillator.set('scattering_length', [ 200.0, 215.0, 225.0, 235.0, 245.0, 255.0, 265.0, 275.0, 285.0, 295.0, 305.0, 315.0, 325.0, 335.0, 345.0, 355.0, 365.0, 375.0, 385.0, 395.0, 405.0, 415.0, 425.0, 435.0, 445.0, 455.0, 465.0, 475.0, 485.0, 495.0, 505.0, 515.0, 525.0, 535.0, 545.0, 555.0, 565.0, 575.0, 585.0, 595.0, 800.0 ], [ 0.290404, 0.387993, 0.465472, 0.553999, 0.654579, 0.768259, 0.89613, 1.03933, 1.19903, 1.37645, 1.57287, 1.78958, 2.02794, 2.28936, 2.57527, 2.88715, 3.22652, 3.59497, 3.9941, 4.42558, 4.89108, 5.3924, 5.93129, 6.50958, 7.12918, 7.79198, 8.49995, 9.2551, 10.0595, 10.9152, 11.8243, 12.7892, 13.8118, 14.8947, 16.04, 17.25, 18.5273, 19.8742, 21.2933, 22.7871, 22.7871 ]) +labppo_scintillator.set('refractive_index', wavelengths=[ 200.0, 213.0, 215.0, 217.0, 219.0, 221.0, 223.0, 225.0, 227.0, 229.0, 231.0, 233.0, 235.0, 237.0, 239.0, 241.0, 243.0, 245.0, 247.0, 249.0, 251.0, 253.0, 255.0, 257.0, 259.0, 261.0, 263.0, 265.0, 267.0, 269.0, 271.0, 273.0, 275.0, 277.0, 279.0, 281.0, 283.0, 285.0, 287.0, 289.0, 291.0, 293.0, 295.0, 297.0, 299.0, 301.0, 303.0, 305.0, 307.0, 309.0, 311.0, 313.0, 315.0, 317.0, 319.0, 321.0, 323.0, 325.0, 327.0, 329.0, 331.0, 333.0, 335.0, 337.0, 339.0, 341.0, 343.0, 345.0, 347.0, 349.0, 351.0, 353.0, 355.0, 357.0, 359.0, 361.0, 363.0, 365.0, 367.0, 369.0, 371.0, 373.0, 375.0, 377.0, 379.0, 381.0, 383.0, 385.0, 387.0, 389.0, 391.0, 393.0, 395.0, 397.0, 399.0, 401.0, 403.0, 405.0, 407.0, 409.0, 411.0, 413.0, 415.0, 417.0, 419.0, 421.0, 423.0, 425.0, 427.0, 429.0, 431.0, 433.0, 435.0, 437.0, 439.0, 441.0, 443.0, 445.0, 447.0, 449.0, 451.0, 453.0, 455.0, 457.0, 459.0, 461.0, 463.0, 465.0, 467.0, 469.0, 471.0, 473.0, 475.0, 477.0, 479.0, 481.0, 483.0, 485.0, 487.0, 489.0, 491.0, 493.0, 495.0, 497.0, 499.0, 501.0, 503.0, 505.0, 507.0, 509.0, 511.0, 513.0, 515.0, 517.0, 519.0, 521.0, 523.0, 525.0, 527.0, 529.0, 531.0, 533.0, 535.0, 537.0, 539.0, 541.0, 543.0, 545.0, 547.0, 549.0, 551.0, 553.0, 555.0, 557.0, 559.0, 561.0, 563.0, 565.0, 567.0, 569.0, 571.0, 573.0, 575.0, 577.0, 579.0, 581.0, 583.0, 585.0, 587.0, 589.0, 591.0, 593.0, 595.0, 597.0, 599.0, 601.0, 603.0, 605.0, 607.0, 609.0, 611.0, 613.0, 615.0, 617.0, 619.0, 621.0, 623.0, 625.0, 627.0, 629.0, 631.0, 633.0, 635.0, 637.0, 639.0, 641.0, 643.0, 645.0, 647.0, 649.0, 651.0, 653.0, 655.0, 657.0, 659.0, 661.0, 663.0, 665.0, 667.0, 669.0, 671.0, 673.0, 675.0, 677.0, 679.0, 681.0, 683.0, 685.0, 687.0, 689.0, 691.0, 693.0, 695.0, 697.0, 699.0, 701.0, 703.0, 705.0, 707.0, 709.0, 711.0, 713.0, 715.0, 717.0, 719.0, 721.0, 723.0, 725.0, 727.0, 729.0, 731.0, 733.0, 735.0, 737.0, 739.0, 741.0, 743.0, 745.0, 747.0, 749.0, 751.0, 753.0, 755.0, 757.0, 759.0, 761.0, 763.0, 765.0, 767.0, 769.0, 771.0, 773.0, 775.0, 777.0, 779.0, 781.0, 783.0, 785.0, 787.0, 789.0, 791.0, 793.0, 795.0, 797.0, 799.0, 800.0 ], value=[ 1.75541, 1.75541, 1.75541, 1.7526, 1.74552, 1.73503, 1.72202, 1.70736, 1.69194, 1.67662, 1.66228, 1.65227, 1.64307, 1.63526, 1.62846, 1.62243, 1.61699, 1.61204, 1.60748, 1.60326, 1.59933, 1.59565, 1.59218, 1.58892, 1.58582, 1.58288, 1.58009, 1.57742, 1.57488, 1.57244, 1.57011, 1.56787, 1.56571, 1.56365, 1.56165, 1.55973, 1.55788, 1.55609, 1.55437, 1.5527, 1.55108, 1.54952, 1.548, 1.54654, 1.54511, 1.54373, 1.54239, 1.54109, 1.53982, 1.53859, 1.5374, 1.53623, 1.5351, 1.534, 1.53293, 1.53188, 1.53086, 1.52987, 1.5289, 1.52795, 1.52703, 1.52613, 1.52525, 1.52439, 1.52355, 1.52273, 1.52193, 1.52115, 1.52039, 1.51964, 1.51891, 1.51819, 1.51749, 1.5168, 1.51613, 1.51547, 1.51483, 1.5142, 1.51358, 1.51297, 1.51238, 1.5118, 1.51122, 1.51066, 1.51012, 1.50958, 1.50905, 1.50853, 1.50802, 1.50752, 1.50703, 1.50655, 1.50608, 1.50561, 1.50516, 1.50471, 1.50427, 1.50383, 1.50341, 1.50299, 1.50258, 1.50218, 1.50178, 1.50139, 1.501, 1.50063, 1.50025, 1.49989, 1.49953, 1.49918, 1.49883, 1.49848, 1.49815, 1.49781, 1.49749, 1.49717, 1.49685, 1.49654, 1.49623, 1.49593, 1.49563, 1.49533, 1.49504, 1.49476, 1.49448, 1.4942, 1.49392, 1.49366, 1.49339, 1.49313, 1.49287, 1.49261, 1.49236, 1.49212, 1.49187, 1.49163, 1.49139, 1.49116, 1.49093, 1.4907, 1.49047, 1.49025, 1.49003, 1.48982, 1.4896, 1.48939, 1.48918, 1.48898, 1.48877, 1.48857, 1.48838, 1.48818, 1.48799, 1.4878, 1.48761, 1.48742, 1.48724, 1.48706, 1.48688, 1.4867, 1.48653, 1.48635, 1.48618, 1.48601, 1.48585, 1.48568, 1.48552, 1.48536, 1.4852, 1.48504, 1.48488, 1.48473, 1.48458, 1.48443, 1.48428, 1.48413, 1.48398, 1.48384, 1.4837, 1.48356, 1.48342, 1.48328, 1.48314, 1.48301, 1.48287, 1.48274, 1.48261, 1.48248, 1.48235, 1.48223, 1.4821, 1.48198, 1.48185, 1.48173, 1.48161, 1.48149, 1.48137, 1.48126, 1.48114, 1.48103, 1.48091, 1.4808, 1.48069, 1.48058, 1.48047, 1.48036, 1.48025, 1.48015, 1.48004, 1.47994, 1.47984, 1.47973, 1.47963, 1.47953, 1.47943, 1.47933, 1.47924, 1.47914, 1.47904, 1.47895, 1.47885, 1.47876, 1.47867, 1.47858, 1.47849, 1.4784, 1.47831, 1.47822, 1.47813, 1.47804, 1.47796, 1.47787, 1.47779, 1.4777, 1.47762, 1.47754, 1.47745, 1.47737, 1.47729, 1.47721, 1.47713, 1.47705, 1.47698, 1.4769, 1.47682, 1.47675, 1.47667, 1.4766, 1.47652, 1.47645, 1.47637, 1.4763, 1.47623, 1.47616, 1.47609, 1.47602, 1.47595, 1.47588, 1.47581, 1.47574, 1.47567, 1.4756, 1.47554, 1.47547, 1.4754, 1.47534, 1.47527, 1.47521, 1.47515, 1.47508, 1.47502, 1.47496, 1.47489, 1.47483, 1.47477, 1.47471, 1.47465, 1.47459, 1.47453, 1.47447, 1.47441, 1.47435, 1.47429, 1.47424, 1.47418, 1.47412, 1.47407, 1.47401, 1.47395, 1.4739, 1.47384, 1.47379, 1.47373, 1.47368, 1.47363, 1.4736 ]) +labppo_scintillator.set('absorption_length', wavelengths=[ 200.0, 215.0, 225.0, 235.0, 245.0, 255.0, 265.0, 275.0, 285.0, 295.0, 305.0, 315.0, 325.0, 335.0, 345.0, 355.0, 365.0, 375.0, 385.0, 395.0, 405.0, 415.0, 425.0, 435.0, 445.0, 455.0, 465.0, 475.0, 485.0, 495.0, 505.0, 515.0, 525.0, 535.0, 545.0, 555.0, 565.0, 575.0, 585.0, 595.0, 800.0 ], value=[ 1.54944e-07, 2.99637e-07, 1.72048e-05, 3.41222e-05, 1.53692e-05, 7.92675e-06, 1.00564e-05, 0.000270688, 0.000949206, 0.00134987, 0.00322232, 0.011032, 0.0379538, 0.708042, 2.68045, 4.20713, 5.71163, 7.92686, 9.90035, 17.239, 28.2318, 41.8939, 67.3786, 120.517, 347.953, 412.636, 147.016, 80.7207, 54.7561, 45.1483, 39.4147, 39.1438, 42.199, 41.6337, 35.783, 47.614, 58.3934, 43.0845, 32.0874, 35.6067, 35.6067 ]) +labppo_scintillator.set('absorption_wls_length', wavelengths=[ 200.0, 215.0, 225.0, 235.0, 245.0, 255.0, 265.0, 275.0, 285.0, 295.0, 305.0, 315.0, 325.0, 335.0, 345.0, 355.0, 365.0, 375.0, 385.0, 395.0, 405.0, 415.0, 425.0, 435.0, 445.0, 455.0, 465.0, 475.0, 485.0, 495.0, 505.0, 515.0, 525.0, 535.0, 545.0, 555.0, 565.0, 575.0, 585.0, 595.0, 800.0 ], value=[ 0.00109808, 0.000693615, 0.000437277, 0.000276096, 0.000173933, 0.000103044, 6.03065e-05, 3.56199e-05, 2.25651e-05, 1.60844e-05, 1.62228e-05, 2.21832e-05, 3.03067e-05, 5.24247e-05, 0.000484492, 0.0127541, 0.294651, 4.15892, 15.6626, 25.4821, 36.1763, 50.723, 63.9962, 77.7386, 89.8611, 103.551, 115.859, 125.442, 146.664, 163.748, 194.986, 194.652, 202.262, 234.033, 262.861, 279.883, 323.991, 299.632, 287.367, 303.42, 303.42 ]) +labppo_scintillator.set('scattering_length', wavelengths=[ 200.0, 215.0, 225.0, 235.0, 245.0, 255.0, 265.0, 275.0, 285.0, 295.0, 305.0, 315.0, 325.0, 335.0, 345.0, 355.0, 365.0, 375.0, 385.0, 395.0, 405.0, 415.0, 425.0, 435.0, 445.0, 455.0, 465.0, 475.0, 485.0, 495.0, 505.0, 515.0, 525.0, 535.0, 545.0, 555.0, 565.0, 575.0, 585.0, 595.0, 800.0 ], value=[ 0.290404, 0.387993, 0.465472, 0.553999, 0.654579, 0.768259, 0.89613, 1.03933, 1.19903, 1.37645, 1.57287, 1.78958, 2.02794, 2.28936, 2.57527, 2.88715, 3.22652, 3.59497, 3.9941, 4.42558, 4.89108, 5.3924, 5.93129, 6.50958, 7.12918, 7.79198, 8.49995, 9.2551, 10.0595, 10.9152, 11.8243, 12.7892, 13.8118, 14.8947, 16.04, 17.25, 18.5273, 19.8742, 21.2933, 22.7871, 22.7871 ]) -- cgit From 142e790aeae55d405f7921527cd8e869d36c5671 Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Wed, 10 Aug 2011 19:22:06 -0400 Subject: Use WCSim properties for water with the lbne detector --- detectors/lbne.py | 6 ++--- optics.py | 74 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ sim.py | 4 +-- 3 files changed, 79 insertions(+), 5 deletions(-) diff --git a/detectors/lbne.py b/detectors/lbne.py index 3e8ae7b..3e27907 100644 --- a/detectors/lbne.py +++ b/detectors/lbne.py @@ -14,16 +14,16 @@ from make import cylinder def build_lbne(radius, height, nstrings, pmts_per_string, endcap_spacing, physical_model=True): if physical_model: - pmt = build_12inch_pmt_with_lc() + pmt = build_12inch_pmt_with_lc(outer_material=water_wcsim) else: - pmt = build_12inch_pmt_shell() + pmt = build_12inch_pmt_shell(outer_material=water_wcsim) lbne = Geometry() # outer cylinder cylinder_mesh = cylinder(radius, radius, height+height/(pmts_per_string-1), theta=(2*np.pi/nstrings)/4) cylinder_mesh.vertices = rotate(cylinder_mesh.vertices, np.pi/2, (-1,0,0)) - lbne.add_solid(Solid(cylinder_mesh, water, vacuum, black_surface, 0xff0000ff)) + lbne.add_solid(Solid(cylinder_mesh, water_wcsim, vacuum, black_surface, 0xff0000ff)) lbne.pmtids = [] diff --git a/optics.py b/optics.py index d65c5f9..c31d209 100644 --- a/optics.py +++ b/optics.py @@ -47,6 +47,8 @@ r7081hqe_photocathode.absorb = r7081hqe_photocathode.detect # remaining photons are diffusely reflected r7081hqe_photocathode.set('reflect_diffuse', 1.0 - r7081hqe_photocathode.detect[:,1] - r7081hqe_photocathode.absorb[:,1], wavelengths=r7081hqe_photocathode.detect[:,0]) +######################## SNO+ materials ############################## + # water data comes from 'lightwater_sno' material in the SNO+ optics database water = Material('water') water.density = 1.0 # g/cm^3 @@ -168,3 +170,75 @@ labppo_scintillator.set('absorption_length', wavelengths=[ 200.0, 215.0, 225.0, labppo_scintillator.set('absorption_wls_length', wavelengths=[ 200.0, 215.0, 225.0, 235.0, 245.0, 255.0, 265.0, 275.0, 285.0, 295.0, 305.0, 315.0, 325.0, 335.0, 345.0, 355.0, 365.0, 375.0, 385.0, 395.0, 405.0, 415.0, 425.0, 435.0, 445.0, 455.0, 465.0, 475.0, 485.0, 495.0, 505.0, 515.0, 525.0, 535.0, 545.0, 555.0, 565.0, 575.0, 585.0, 595.0, 800.0 ], value=[ 0.00109808, 0.000693615, 0.000437277, 0.000276096, 0.000173933, 0.000103044, 6.03065e-05, 3.56199e-05, 2.25651e-05, 1.60844e-05, 1.62228e-05, 2.21832e-05, 3.03067e-05, 5.24247e-05, 0.000484492, 0.0127541, 0.294651, 4.15892, 15.6626, 25.4821, 36.1763, 50.723, 63.9962, 77.7386, 89.8611, 103.551, 115.859, 125.442, 146.664, 163.748, 194.986, 194.652, 202.262, 234.033, 262.861, 279.883, 323.991, 299.632, 287.367, 303.42, 303.42 ]) labppo_scintillator.set('scattering_length', wavelengths=[ 200.0, 215.0, 225.0, 235.0, 245.0, 255.0, 265.0, 275.0, 285.0, 295.0, 305.0, 315.0, 325.0, 335.0, 345.0, 355.0, 365.0, 375.0, 385.0, 395.0, 405.0, 415.0, 425.0, 435.0, 445.0, 455.0, 465.0, 475.0, 485.0, 495.0, 505.0, 515.0, 525.0, 535.0, 545.0, 555.0, 565.0, 575.0, 585.0, 595.0, 800.0 ], value=[ 0.290404, 0.387993, 0.465472, 0.553999, 0.654579, 0.768259, 0.89613, 1.03933, 1.19903, 1.37645, 1.57287, 1.78958, 2.02794, 2.28936, 2.57527, 2.88715, 3.22652, 3.59497, 3.9941, 4.42558, 4.89108, 5.3924, 5.93129, 6.50958, 7.12918, 7.79198, 8.49995, 9.2551, 10.0595, 10.9152, 11.8243, 12.7892, 13.8118, 14.8947, 16.04, 17.25, 18.5273, 19.8742, 21.2933, 22.7871, 22.7871 ]) +################### WCSim materials ##################### + +water_wcsim = Material('water') +water_wcsim.density = 1.0 # g/cm^3 +water_wcsim.composition = { 'H' : 0.1119, 'O' : 0.8881 } # fraction by mass +hc_over_GeV = 1.2398424468024265e-06 # h_Planck * c_light / GeV / nanometer +wcsim_wavelengths = hc_over_GeV / np.array([ 1.56962e-09, 1.58974e-09, 1.61039e-09, 1.63157e-09, + 1.65333e-09, 1.67567e-09, 1.69863e-09, 1.72222e-09, + 1.74647e-09, 1.77142e-09,1.7971e-09, 1.82352e-09, + 1.85074e-09, 1.87878e-09, 1.90769e-09, 1.93749e-09, + 1.96825e-09, 1.99999e-09, 2.03278e-09, 2.06666e-09, + 2.10169e-09, 2.13793e-09, 2.17543e-09, 2.21428e-09, + 2.25454e-09, 2.29629e-09, 2.33962e-09, 2.38461e-09, + 2.43137e-09, 2.47999e-09, 2.53061e-09, 2.58333e-09, + 2.63829e-09, 2.69565e-09, 2.75555e-09, 2.81817e-09, + 2.88371e-09, 2.95237e-09, 3.02438e-09, 3.09999e-09, + 3.17948e-09, 3.26315e-09, 3.35134e-09, 3.44444e-09, + 3.54285e-09, 3.64705e-09, 3.75757e-09, 3.87499e-09, + 3.99999e-09, 4.13332e-09, 4.27585e-09, 4.42856e-09, + 4.59258e-09, 4.76922e-09, 4.95999e-09, 5.16665e-09, + 5.39129e-09, 5.63635e-09, 5.90475e-09, 6.19998e-09 ])[::-1] #reversed + +water_wcsim.set('refractive_index', + wavelengths=wcsim_wavelengths, + value=np.array([1.32885, 1.32906, 1.32927, 1.32948, 1.3297, 1.32992, 1.33014, + 1.33037, 1.3306, 1.33084, 1.33109, 1.33134, 1.3316, 1.33186, 1.33213, + 1.33241, 1.3327, 1.33299, 1.33329, 1.33361, 1.33393, 1.33427, 1.33462, + 1.33498, 1.33536, 1.33576, 1.33617, 1.3366, 1.33705, 1.33753, 1.33803, + 1.33855, 1.33911, 1.3397, 1.34033, 1.341, 1.34172, 1.34248, 1.34331, + 1.34419, 1.34515, 1.3462, 1.34733, 1.34858, 1.34994, 1.35145, 1.35312, + 1.35498, 1.35707, 1.35943, 1.36211, 1.36518, 1.36872, 1.37287, 1.37776, + 1.38362, 1.39074, 1.39956, 1.41075, 1.42535])[::-1] #reversed +) +water_wcsim.set('absorption_length', + wavelengths=wcsim_wavelengths, + value=np.array([22.8154, 28.6144, 35.9923, 45.4086, 57.4650, + 72.9526, 75, 81.2317, 120.901, 160.243, + 193.797, 215.045, 227.786, 243.893, 294.113, + 321.735, 342.931, 362.967, 378.212, 449.602, + 740.143, 1116.06, 1438.78, 1615.48, 1769.86, + 2109.67, 2304.13, 2444.97, 3076.83, 4901.5, + 6666.57, 7873.95, 9433.81, 10214.5, 10845.8, + 15746.9, 20201.8, 22025.8, 21142.2, 15083.9, + 11751, 8795.34, 8741.23, 7102.37, 6060.68, + 4498.56, 3039.56, 2232.2, 1938, 1811.58, + 1610.32, 1338.7, 1095.3, 977.525, 965.258, + 1082.86, 876.434, 633.723, 389.87, 142.011])[::-1] / 100.0 # reversed, cm->m + ) + +water_wcsim.set('scattering_length', + wavelengths=wcsim_wavelengths, + value=np.array([167024.4, 158726.7, 150742, + 143062.5, 135680.2, 128587.4, + 121776.3, 115239.5, 108969.5, + 102958.8, 97200.35, 91686.86, + 86411.33, 81366.79, 76546.42, + 71943.46, 67551.29, 63363.36, + 59373.25, 55574.61, 51961.24, + 48527.00, 45265.87, 42171.94, + 39239.39, 36462.50, 33835.68, + 31353.41, 29010.30, 26801.03, + 24720.42, 22763.36, 20924.88, + 19200.07, 17584.16, 16072.45, + 14660.38, 13343.46, 12117.33, + 10977.70, 9920.416, 8941.407, + 8036.711, 7202.470, 6434.927, + 5730.429, 5085.425, 4496.467, + 3960.210, 3473.413, 3032.937, + 2635.746, 2278.907, 1959.588, + 1675.064, 1422.710, 1200.004, + 1004.528, 833.9666, 686.1063])[::-1] / 100.0 * 0.625 # reversed, cm -> m, * magic tuning constant + ) diff --git a/sim.py b/sim.py index 7644266..9f89cd1 100755 --- a/sim.py +++ b/sim.py @@ -120,14 +120,14 @@ def main(): detector.build(bits=options.nbits) print >>sys.stderr, 'Creating generator...' - detector_material = optics.water + detector_material = optics.water_wcsim generator_thread = GeneratorProcess(particle=options.particle, energy=options.energy, position=position, direction=direction, nevents=options.nevents, material=detector_material) - print >>sys.stderr, 'WARNING: ASSUMING DETECTOR IS WATER!!' + print >>sys.stderr, 'WARNING: ASSUMING DETECTOR IS WCSIM WATER!!' # Do this now so we can get ahead of the photon propagation print >>sys.stderr, 'Starting GEANT4 generator...' -- cgit From 56afb978b2416ee9a14ecacdf41ab996d1747b66 Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Wed, 10 Aug 2011 20:56:14 -0400 Subject: Set the GEANT4 and CUDA RNG seeds using current time and process ID if not set on command line. --- g4gen.py | 7 ++++++- sim.py | 26 +++++++++++++++++++------- 2 files changed, 25 insertions(+), 8 deletions(-) diff --git a/g4gen.py b/g4gen.py index 409570f..718df6e 100644 --- a/g4gen.py +++ b/g4gen.py @@ -18,7 +18,7 @@ except: import G4chroma class G4Generator(object): - def __init__(self, material): + def __init__(self, material, seed=None): '''Create generator to produce photons inside the specified material. material: chroma.geometry.Material object with density, @@ -26,7 +26,12 @@ class G4Generator(object): composition dictionary should be { element_symbol : fraction_by_weight, ... }. + seed: Random number generator seed for HepRandom. If None, + generator is not seeded. ''' + if seed is not None: + HepRandom.setTheSeed(seed) + g4py.NISTmaterials.Construct() g4py.ezgeom.Construct() self.physics_list = G4chroma.ChromaPhysicsList() diff --git a/sim.py b/sim.py index 9f89cd1..26a6b68 100755 --- a/sim.py +++ b/sim.py @@ -2,6 +2,7 @@ import sys import optparse import time +import os import multiprocessing import detectors @@ -13,6 +14,11 @@ import numpy as np import math import ROOT +def pick_seed(): + '''Returns a seed for a random number generator selected using + a mixture of the current time and the current process ID.''' + return int(time.time()) ^ (os.getpid() << 16) + def info(type, value, tb): if hasattr(sys, 'ps1') or not sys.stderr.isatty(): # we are in interactive mode or we don't have a tty-like @@ -26,10 +32,8 @@ def info(type, value, tb): # ...then start the debugger in post-mortem mode. pdb.pm() -sys.excepthook = info - class GeneratorProcess(multiprocessing.Process): - def __init__(self, particle, energy, position, direction, nevents, material): + def __init__(self, particle, energy, position, direction, nevents, material, seed=None): multiprocessing.Process.__init__(self) self.particle = particle @@ -38,11 +42,12 @@ class GeneratorProcess(multiprocessing.Process): self.direction = direction self.nevents = nevents self.material = material + self.seed = seed self.queue = multiprocessing.Queue() def run(self): print >>sys.stderr, 'Starting generator thread...' - generator = g4gen.G4Generator(self.material) + generator = g4gen.G4Generator(self.material, seed=self.seed) for i in xrange(self.nevents): photons = generator.generate_photons(particle_name=self.particle, @@ -87,7 +92,8 @@ def main(): parser.add_option('-b', type='int', dest='nbits', default=10) parser.add_option('-j', type='int', dest='device', default=None) parser.add_option('-n', type='int', dest='nblocks', default=64) - + parser.add_option('-s', type='int', dest='seed', default=None, + help='Set random number generator seed') parser.add_option('--detector', type='string', dest='detector', default='microlbne') parser.add_option('--nevents', type='int', dest='nevents', default=100) parser.add_option('--particle', type='string', dest='particle', default='e-') @@ -115,6 +121,10 @@ def main(): position = np.array(eval(options.pos), dtype=float) direction = np.array(eval(options.dir), dtype=float) detector = detectors.find(options.detector) + if options.seed is None: + options.seed = pick_seed() + + print >>sys.stderr, 'RNG seed:', options.seed print >>sys.stderr, 'Creating BVH for detector "%s" with %d bits...' % (options.detector, options.nbits) detector.build(bits=options.nbits) @@ -126,7 +136,8 @@ def main(): position=position, direction=direction, nevents=options.nevents, - material=detector_material) + material=detector_material, + seed=options.seed) print >>sys.stderr, 'WARNING: ASSUMING DETECTOR IS WCSIM WATER!!' # Do this now so we can get ahead of the photon propagation @@ -140,7 +151,7 @@ def main(): gpu_worker.load_geometry(detector) print >>sys.stderr, 'Initializing random numbers generators...' - gpu_worker.setup_propagate() + gpu_worker.setup_propagate(seed=options.seed) gpu_worker.setup_daq(max(detector.pmtids)) # Create output file @@ -191,4 +202,5 @@ def main(): print >>sys.stderr, 'Done. %1.1f events/sec, %1.0f photons/sec.' % (options.nevents/(end_sim - start_sim), nphotons/(end_sim - start_sim)) if __name__ == '__main__': + sys.excepthook = info main() -- cgit From 6310592deb66d75f473c63fe286187ca41dae0f6 Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Thu, 11 Aug 2011 11:49:54 -0400 Subject: Switch from texture to float3 array for upper and lower bounds. 10% speed boost! --- gpu.py | 30 ++++++++++++------------------ src/mesh.h | 12 +++++++----- 2 files changed, 19 insertions(+), 23 deletions(-) diff --git a/gpu.py b/gpu.py index e6b6856..ffa17a6 100644 --- a/gpu.py +++ b/gpu.py @@ -143,37 +143,31 @@ class GPU(object): triangles['w'] = ((geometry.material1_index & 0xff) << 24) | ((geometry.material2_index & 0xff) << 16) | ((geometry.surface_index & 0xff) << 8) self.triangles_gpu = gpuarray.to_gpu(triangles) - lower_bounds_float4 = np.empty(geometry.lower_bounds.shape[0], dtype=gpuarray.vec.float4) - lower_bounds_float4['x'] = geometry.lower_bounds[:,0] - lower_bounds_float4['y'] = geometry.lower_bounds[:,1] - lower_bounds_float4['z'] = geometry.lower_bounds[:,2] - self.lower_bounds_gpu = gpuarray.to_gpu(lower_bounds_float4) - - upper_bounds_float4 = np.empty(geometry.upper_bounds.shape[0], dtype=gpuarray.vec.float4) - upper_bounds_float4['x'] = geometry.upper_bounds[:,0] - upper_bounds_float4['y'] = geometry.upper_bounds[:,1] - upper_bounds_float4['z'] = geometry.upper_bounds[:,2] - self.upper_bounds_gpu = gpuarray.to_gpu(upper_bounds_float4) + lower_bounds_float3 = np.empty(geometry.lower_bounds.shape[0], dtype=gpuarray.vec.float3) + lower_bounds_float3['x'] = geometry.lower_bounds[:,0] + lower_bounds_float3['y'] = geometry.lower_bounds[:,1] + lower_bounds_float3['z'] = geometry.lower_bounds[:,2] + self.lower_bounds_gpu = gpuarray.to_gpu(lower_bounds_float3) + + upper_bounds_float3 = np.empty(geometry.upper_bounds.shape[0], dtype=gpuarray.vec.float3) + upper_bounds_float3['x'] = geometry.upper_bounds[:,0] + upper_bounds_float3['y'] = geometry.upper_bounds[:,1] + upper_bounds_float3['z'] = geometry.upper_bounds[:,2] + self.upper_bounds_gpu = gpuarray.to_gpu(upper_bounds_float3) self.colors_gpu = gpuarray.to_gpu(geometry.colors.astype(np.uint32)) self.node_map_gpu = gpuarray.to_gpu(geometry.node_map.astype(np.uint32)) self.node_map_end_gpu = gpuarray.to_gpu(geometry.node_map_end.astype(np.uint32)) self.solid_id_map_gpu = gpuarray.to_gpu(geometry.solid_id.astype(np.uint32)) - self.geo_funcs.set_global_mesh_variables(self.triangles_gpu, self.vertices_gpu, self.colors_gpu, np.uint32(geometry.node_map.size-1), np.uint32(geometry.first_node), block=(1,1,1), grid=(1,1)) + self.geo_funcs.set_global_mesh_variables(self.triangles_gpu, self.vertices_gpu, self.colors_gpu, np.uint32(geometry.node_map.size-1), np.uint32(geometry.first_node), self.lower_bounds_gpu, self.upper_bounds_gpu, block=(1,1,1), grid=(1,1)) - self.lower_bounds_tex = self.module.get_texref('lower_bounds') - self.upper_bounds_tex = self.module.get_texref('upper_bounds') self.node_map_tex = self.module.get_texref('node_map') self.node_map_end_tex = self.module.get_texref('node_map_end') - self.lower_bounds_tex.set_address(self.lower_bounds_gpu.gpudata, self.lower_bounds_gpu.nbytes) - self.upper_bounds_tex.set_address(self.upper_bounds_gpu.gpudata, self.upper_bounds_gpu.nbytes) self.node_map_tex.set_address(self.node_map_gpu.gpudata, self.node_map_gpu.nbytes) self.node_map_end_tex.set_address(self.node_map_end_gpu.gpudata, self.node_map_end_gpu.nbytes) - self.lower_bounds_tex.set_format(cuda.array_format.FLOAT, 4) - self.upper_bounds_tex.set_format(cuda.array_format.FLOAT, 4) self.node_map_tex.set_format(cuda.array_format.UNSIGNED_INT32, 1) self.node_map_end_tex.set_format(cuda.array_format.UNSIGNED_INT32, 1) diff --git a/src/mesh.h b/src/mesh.h index b4714c4..bb30bef 100644 --- a/src/mesh.h +++ b/src/mesh.h @@ -13,8 +13,8 @@ __device__ unsigned int g_start_node; __device__ unsigned int g_first_node; /* lower/upper bounds for the bounding box associated with each node/leaf */ -texture upper_bounds; -texture lower_bounds; +__device__ float3 *g_lower_bounds; +__device__ float3 *g_upper_bounds; /* map to child node/triangle indices */ texture node_map; @@ -38,8 +38,8 @@ __device__ int convert(int c) intersects the bounding box return true, else return false. */ __device__ bool intersect_node(const float3 &origin, const float3 &direction, const int &i) { - float3 lower_bound = make_float3(tex1Dfetch(lower_bounds, i)); - float3 upper_bound = make_float3(tex1Dfetch(upper_bounds, i)); + float3 lower_bound = g_lower_bounds[i]; + float3 upper_bound = g_upper_bounds[i]; return intersect_box(origin, direction, lower_bound, upper_bound); } @@ -134,13 +134,15 @@ __device__ int intersect_mesh(const float3 &origin, const float3& direction, flo extern "C" { -__global__ void set_global_mesh_variables(uint4 *triangles, float3 *vertices, unsigned int *colors, unsigned int start_node, unsigned int first_node) + __global__ void set_global_mesh_variables(uint4 *triangles, float3 *vertices, unsigned int *colors, unsigned int start_node, unsigned int first_node, float3 *lower_bounds, float3 *upper_bounds) { g_triangles = triangles; g_vertices = vertices; g_colors = colors; g_start_node = start_node; g_first_node = first_node; + g_lower_bounds = lower_bounds; + g_upper_bounds = upper_bounds; } __global__ void set_colors(unsigned int *colors) -- cgit From f184de9a5ab72188e03149a7d63a7d0da00a868e Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Thu, 11 Aug 2011 14:20:59 -0400 Subject: Show number of registers used in CUDA kernels --- gpu.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gpu.py b/gpu.py index ffa17a6..16e6f4a 100644 --- a/gpu.py +++ b/gpu.py @@ -72,7 +72,7 @@ class GPU(object): device = cuda.Device(device_id) self.context = device.make_context() print 'device %s' % self.context.get_device().name() - cuda_options = ['-I' + src.dir, '--use_fast_math'] + cuda_options = ['-I' + src.dir, '--use_fast_math', '--ptxas-options=-v'] self.module = SourceModule(src.kernel, options=cuda_options, no_extern_c=True) self.geo_funcs = CUDAFuncs(self.module, ['set_wavelength_range', 'set_material', 'set_surface', 'set_global_mesh_variables', 'color_solids']) -- cgit From 1b9aef380e629c136e045412621483ebb4a4164f Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Thu, 11 Aug 2011 14:21:25 -0400 Subject: No need for __noinline__ now that kernel caching works. --- src/kernel.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kernel.cu b/src/kernel.cu index 17b829c..fe518f6 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -17,7 +17,7 @@ __device__ void fAtomicAdd(float *addr, float data) } } -__device__ __noinline__ void to_diffuse(Photon &p, State &s, curandState &rng, const int &max_steps) +__device__ void to_diffuse(Photon &p, State &s, curandState &rng, const int &max_steps) { int steps = 0; while (steps < max_steps) -- cgit