Commits

Stan Seibert committed 10bad4e

Initial commit of photon detector geometry for LBNE with 5 kiloton (fiducial)
cryostat.

  • Participants

Comments (0)

Files changed (8)

+syntax: glob
+*.pyc
+Chroma_PDS.egg-info

File chroma_pds/__init__.py

Empty file added.

File chroma_pds/detector.py

+'''Chroma geometry describing photon detector system in LBNE.'''
+
+import numpy as np
+from chroma_pds.optics import vacuum, acrylic, liquid_Ar, tpb, \
+                              photocathode, steel, aluminum
+from chroma.geometry import Solid
+from chroma.detector import Detector
+from chroma.make import box
+
+
+#### Basic detector properties
+APA_HEIGHT = 7000.0  # mm
+
+APA_WIDTH = 2520.0  # mm
+APA_SPACING = 4640.0  # mm
+LONGITUDINAL_APA_GAP = 20.0  # mm
+VERTICAL_APA_GAP = 25.0  # mm
+
+INNER_HEIGHT = 16000.0  # mm
+INNER_WIDTH = 15600.0  # mm
+INNER_LENGTH = 27000.0  # mm
+####
+
+
+def build_cryostat(ev_display=False):
+    '''Build the rectangular steel cryostat.'''
+    tank_mesh = box(INNER_LENGTH, INNER_WIDTH, INNER_HEIGHT)
+    if ev_display:
+        color = 0xFF00FFAA
+    else:
+        color = 0xAA00FFAA
+    tank = Solid(tank_mesh, liquid_Ar, vacuum, surface=steel,
+                 color=color)
+
+    return tank
+
+
+def build_paddle(ev_display=False, mask_side=None, double_side=False):
+    '''Build an acrylic bar with optional mask and photo-sensitive detectors
+    on each end.
+
+    Parameters:
+        ev_display: If true, make bars partially transparent for display.
+        mask_side: If None, no mask.  Otherwise -1 or +1 to select the
+                   side for the aluminum mask.
+        double_side: If True, put photo-sensitive detector on both ends
+                     of bar, otherwise just on one end.
+    '''
+    length = 2000.0 + 250.0  # mm (including 90 twist length here)
+    width = 88.0  # mm
+    thickness = 4.8  # mm
+
+    mesh = box(length, thickness, width)
+
+    if ev_display:
+        tpb_color = 0xEEEEEEEE
+        reflector_color = 0xEEEE00EE
+        detector_color = tpb_color
+    else:
+        tpb_color = 0x44EE0000
+        reflector_color = 0x44EEEEEE
+        detector_color = 0x44FFFF00
+
+    surface = np.empty(shape=len(mesh.triangles), dtype=np.object)
+    color = np.empty(shape=len(mesh.triangles), dtype=np.uint32)
+
+    # Pick out the triangles on one end of the paddle for photon detection
+    triangles = mesh.assemble(group=True)
+    triangle_u = triangles[:, 1, :] - triangles[:, 0, :]
+    triangle_v = triangles[:, 2, :] - triangles[:, 0, :]
+    triangle_normal = np.cross(triangle_u, triangle_v)
+    if double_side:
+        pc_side = np.abs(triangle_normal[:, 0]) > 0.01
+    else:
+        pc_side = triangle_normal[:, 0] > 0.01
+
+    surface[~pc_side] = tpb
+    surface[pc_side] = photocathode
+    color[~pc_side] = tpb_color
+    color[pc_side] = detector_color
+
+    # Add a mask if requested
+    if mask_side is not None:
+        assert mask_side in (-1, 1)
+        mask_triangles = mask_side * triangle_normal[:, 1] > 0.01
+        surface[mask_triangles] = aluminum
+        color[mask_triangles] = reflector_color
+
+    paddle = Solid(mesh, acrylic, liquid_Ar, surface=surface,
+                   color=color)
+    return paddle
+
+
+def add_apa(geo, displacement, ev_display=False, num_paddles=10,
+            masked_bars=False, double_side=False):
+    '''Add an anode plane assembly to geo with the given displacement.
+
+    Parameters:
+        geo: chroma.geometry.Geometry
+        displacement: numpy.ndarray(shape=3)
+        ev_display: If true, make bars very translucent for ev_display
+        num_paddles: number of acrylic bars per APA
+        masked_bars: Cover one side of each bar with aluminum.  Alternates on each bar.
+        double_side: Put photo-sensitive detector on each end of bar.
+    '''
+    height = APA_HEIGHT
+    width = APA_WIDTH
+    bar_width = 90.0
+    stride = height / num_paddles
+
+    if ev_display:
+        bar_color = 0xE8EEEEEE
+    else:
+        bar_color = 0x00888888
+
+    #Steel frame components
+    left_bar = Solid(box(bar_width, bar_width, height),
+                     liquid_Ar, liquid_Ar,
+                     surface=steel, color=bar_color)
+    left_bar.mesh.vertices[:, 0] -= (width - bar_width) / 2
+
+    right_bar = Solid(box(bar_width, bar_width, height),
+                      liquid_Ar, liquid_Ar,
+                      surface=steel, color=bar_color)
+    right_bar.mesh.vertices[:, 0] += (width - bar_width) / 2
+
+    top_bar = Solid(box(width - bar_width, bar_width, bar_width),
+                    liquid_Ar, liquid_Ar,
+                    surface=steel, color=bar_color)
+    top_bar.mesh.vertices[:, 2] -= (height - bar_width) / 2
+
+    bottom_bar = Solid(box(width - bar_width, bar_width, bar_width),
+                       liquid_Ar, liquid_Ar,
+                       surface=steel, color=bar_color)
+    bottom_bar.mesh.vertices[:, 2] += (height - bar_width) / 2
+
+    # Assemble APA frame
+    apa = left_bar + right_bar + top_bar + bottom_bar
+    geo.add_solid(apa, displacement=displacement)
+
+    # Add acrylic bars
+    for i, z in enumerate(np.arange(-(height - stride) / 2,
+                                    (height - stride) / 2 + 0.001,
+                                    stride)):
+        if masked_bars:
+            mask_side = (i % 2) * 2 - 1  # alternate -1 and 1
+        else:
+            mask_side = None
+        paddle = build_paddle(ev_display=ev_display, mask_side=mask_side,
+                              double_side=double_side)
+        paddle.mesh.vertices[:] += (-bar_width / 4, 0, z)
+        geo.add_pmt(paddle, displacement=displacement)
+
+
+def make_cathode_plane(surface, y, height, length, ev_display=False):
+    '''Create a solid cathode plane (infinitesmial thickness).
+
+    Parameters:
+        surface: chroma.geometry.Surface
+        y: Offset of plane in y direction (mm)
+        height: Height of plane (mm)
+        ev_display: If True, make plane render transparent for display
+    '''
+    center = (0, y, 0)
+    thickness = 3.0  # mm
+    plate = box(length, thickness, height, center)
+    if ev_display:
+        color = 0xBB222244
+    else:
+        color = 0x44222244
+    return Solid(plate, liquid_Ar, liquid_Ar, surface=surface, color=color)
+
+
+def lar5(ev_display=False, efficiency_scale=1.0, masked_bars=False,
+         double_side=False):
+    '''Create an LBNE cryostat (5 kton fiducial).
+
+    Parameters:
+        ev_display: If True, sets the colors of most detector elements
+            to transparent for event display.
+        efficiency_scale: Multiply the photocathode efficiency by a factor
+        masked_bars: If True, make acrylic bars have one side masked, alternating
+            sides down the APA.
+        double_side: If True, put photo-sensitive detectors on both ends of bars.
+
+    Returns: chroma.detector.Detector
+    '''
+    photocathode.detect[:, 1] *= efficiency_scale
+    geo = Detector(liquid_Ar)
+
+    geo.add_solid(build_cryostat(ev_display=ev_display))
+    nplanes = 3
+    nlong = 10
+    apa_count = 0
+
+    # Anode planes
+    for x in np.arange(-(APA_WIDTH + LONGITUDINAL_APA_GAP) * (nlong - 1) / 2.0,
+                       (APA_WIDTH + LONGITUDINAL_APA_GAP) * (nlong - 1) / 2.0 + 0.001,
+                        APA_WIDTH + LONGITUDINAL_APA_GAP):
+        for y in np.arange(-(APA_SPACING) * (nplanes - 1) / 2.0,
+                            (APA_SPACING) * (nplanes - 1) / 2.0 + 0.001,
+                            APA_SPACING):
+            for z in [-(APA_HEIGHT + VERTICAL_APA_GAP) / 2,
+                       (APA_HEIGHT + VERTICAL_APA_GAP) / 2]:
+                displacement = (x + 500, y, z)
+                add_apa(geo, displacement=displacement, ev_display=ev_display,
+                        masked_bars=masked_bars, double_side=double_side)
+                apa_count += 1
+
+    # Cathode planes
+    for y in np.arange(-(APA_SPACING) * (nplanes) / 2.0,
+                        (APA_SPACING) * (nplanes) / 2.0 + 0.001,
+                        APA_SPACING):
+        geo.add_solid(make_cathode_plane(steel, y,
+                                height=(APA_HEIGHT + VERTICAL_APA_GAP) * 2,
+                                length=(APA_WIDTH + LONGITUDINAL_APA_GAP) * 10))
+
+    geo.set_time_dist_gaussian(0.34, -10, 10)
+    geo.set_charge_dist_gaussian(1.0, 0.01, 0.5, 1.5)
+    print 'APAs:', apa_count
+    return geo
+
+
+def lar5_mask():
+    '''Create a 5 kton cryostat with masked bars.'''
+    return lar5(masked_bars=True)
+
+
+def lar5_mask_double():
+    '''Create a 5 kton cryostat with double-ended photon readout.'''
+    return lar5(masked_bars=True, double_side=True)
+
+
+def lar5_display():
+    '''Create a 5 kton cryostat with triangles colored for event display.'''
+    return lar5(ev_display=True)

File chroma_pds/optics.py

+from chroma.geometry import Surface, Material
+import numpy as np
+
+# For properties that are segmented between EUV and visible
+# with no wavelength-dependent measurements
+basic_wavelengths = [60.0, 200.0, 210.0, 800.0]
+
+
+# Vacuum
+vacuum = Material('vacuum')
+vacuum.set('refractive_index', 1.0)
+vacuum.set('absorption_length', 1e6)
+vacuum.set('scattering_length', 1e6)
+
+# liquid argon
+liquid_Ar = Material('liquid_Ar')
+# 86K data from Phys. Rev. 181, 1297, extrapolated down to 130 nm and up to 800. The 60 nm value is bogus
+liquid_Ar.set('refractive_index',
+        np.array([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], dtype=np.float32),
+        wavelengths=np.array([60, 130, 140, 150, 160, 170, 180, 190, 200, 210,
+            220, 230, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 340,
+            350, 360, 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], dtype=np.float32))
+# extrap from 900 mm at \lambda = 128 nm
+liquid_Ar.set('scattering_length',
+        np.array([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], dtype=np.float32),
+        wavelengths=np.array([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], dtype=np.float32))
+# Guess.  Depends on actual impurities in liquid argon
+liquid_Ar.set('absorption_length', [3000.0, 3000.0, 10000.0, 10000.0],
+              wavelengths=basic_wavelengths)
+
+# Acrylic
+acrylic = Material('acrylic')
+# Needs to be replaced with measurements of LBNE acrylic!
+acrylic.set('refractive_index',
+        np.array([1.597, 1.597, 1.584, 1.573, 1.564, 1.556, 1.550,
+        1.544, 1.539, 1.534, 1.531, 1.527, 1.524, 1.521, 1.519, 1.516, 1.514,
+        1.512, 1.510, 1.509, 1.507, 1.506, 1.505, 1.503, 1.502, 1.501, 1.500,
+        1.499, 1.499, 1.498, 1.497, 1.496, 1.496, 1.495, 1.494, 1.494, 1.493,
+        1.493, 1.492, 1.492, 1.491, 1.491, 1.490, 1.490, 1.490, 1.489, 1.489,
+        1.488, 1.488, 1.488, 1.488, 1.487, 1.487, 1.487, 1.486, 1.486, 1.486,
+        1.486, 1.485, 1.485, 1.485, 1.485], dtype=np.float32),
+        wavelengths=np.array([60, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300,
+        310, 320, 330, 340, 350, 360, 370, 380, 390, 400, 410, 420, 430, 440,
+        450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580,
+        590, 600, 610, 620, 630, 640, 650, 660, 670, 680, 690, 700, 710, 720,
+        730, 740, 750, 760, 770, 780, 790, 800], dtype=np.float32))
+acrylic.set('absorption_length', [0.001, 0.001, 2500.0, 2500.0],
+    wavelengths=[60.0, 380.0, 390.0, 800.0])
+# What is real rayleigh scattering length of acrylic?? (set to infinity for now)
+acrylic.set('scattering_length', [1e9, 1e9, 1e9, 1e9],
+    wavelengths=[60.0, 380.0, 390.0, 800.0])
+
+# TPB (using chroma surface model #2 for reemission)
+tpb = Surface('tpb', model=2)
+tpb.set('reflect_specular', [0.0, 0.0, 0.0, 0.0], wavelengths=basic_wavelengths)
+tpb.set('reflect_diffuse',  [0.0, 0.0, 0.0, 0.0], wavelengths=basic_wavelengths)
+tpb.set('absorb',  [1.0, 1.0, 0.0, 0.0], wavelengths=basic_wavelengths)
+tpb.set('reemit',  [1.0, 1.0, 0.0, 0.0], wavelengths=basic_wavelengths)
+reemisson_wavelengths = np.array([60, 300, 350.0, 351.5, 353.1,
+    354.6, 356.1, 357.7, 359.2, 360.8, 362.3, 363.8, 365.4,
+    366.9, 368.5, 370.0, 371.5, 373.1, 374.6, 376.1, 377.7,
+    379.2, 380.8, 382.3, 383.8, 385.4, 386.9, 388.4, 390.0,
+    391.5, 393.1, 394.6, 396.1, 397.7, 399.2, 400.7, 402.3,
+    403.8, 405.4, 406.9, 408.4, 410.0, 411.5, 413.0, 414.6,
+    416.1, 417.7, 419.2, 420.7, 422.3, 423.8, 425.4, 426.9,
+    428.4, 430.0, 431.5, 433.0, 434.6, 436.1, 437.7, 439.2,
+    440.7, 442.3, 443.8, 445.3, 446.9, 448.4, 450.0, 451.5,
+    453.0, 454.6, 456.1, 457.6, 459.2, 460.7, 462.3, 463.8,
+    465.3, 466.9, 468.4, 469.9, 471.5, 473.0, 474.6, 476.1,
+    477.6, 479.2, 480.7, 482.2, 483.8, 485.3, 486.9, 488.4,
+    489.9, 491.5, 493.0, 494.6, 496.1, 497.6, 499.2, 500.7,
+    502.2, 503.8, 505.3, 506.9, 508.4, 509.9, 511.5, 513.0,
+    514.5, 516.1, 517.6, 519.2, 520.7, 522.2, 523.8, 525.3,
+    526.8, 528.4, 529.9, 531.5, 533.0, 534.5, 536.1, 537.6,
+    539.1, 540.7, 542.2, 543.8, 545.3, 546.8, 548.4, 549.9,
+    551.4, 553.0, 554.5, 556.1, 557.6, 559.1, 560.7, 562.2,
+    563.8, 565.3, 566.8, 568.4, 569.9, 571.4, 573.0, 574.5,
+    576.1, 577.6, 579.1, 580.7, 582.2, 583.7, 585.3, 586.8,
+    588.4, 589.9, 591.4, 593.0, 594.5, 596.0, 597.6, 599.1], dtype=np.float32)
+reemission_spec = np.array([0, 0, 0.0, 1.0, 0.0, 0.9,
+    0.0, 0.7, 0.2, 0.8, 0.1, 0.3, 0.0, 0.8, 0.4,
+    3.7, 1.2, 1.3, 0.0, 0.5, 0.0, 0.8, 0.0, 1.8,
+    0.1, 2.4, 1.9, 0.7, 2.6, 3.8, 6.5, 6.4, 9.9,
+    10.3, 13.8, 14.1, 17.3, 18.5, 19.6, 20.8, 22.2,
+    22.8, 23.0, 25.3, 26.5, 26.5, 26.2, 27.4, 27.2,
+    27.0, 26.5, 25.5, 25.3, 25.4, 25.1, 23.6, 24.5,
+    23.0, 23.5, 22.9, 22.4, 22.3, 22.0, 21.8, 21.0,
+    20.7, 19.5, 19.7, 17.4, 18.5, 16.0, 16.9, 15.1,
+    16.3, 14.9, 15.2, 16.5, 13.8, 11.7, 10.8, 10.6,
+    11.7, 9.3, 7.5, 9.6, 9.1, 7.6, 7.7, 7.6, 7.7,
+    7.5, 7.5, 6.2, 5.4, 5.5, 6.2, 5.7, 3.6, 4.1,
+    4.8, 4.6, 5.8, 4.8, 4.0, 7.0, 7.6, 3.5, 2.8,
+    3.4, 4.4, 2.9, 3.6, 2.5, 4.3, 2.0, 1.9, 2.0,
+    2.0, 1.5, 2.4, 2.6, 0.0, 3.5, 0.0, 0.8, 0.0,
+    0.8, 1.8, 0.9, 3.1, 2.6, 2.1, 0.8, 1.2, 1.7,
+    1.0, 0.9, 0.5, 1.4, 1.1, 1.2, 1.0, 0.7, 0.4,
+    1.3, 0.3, 0.0, 0.7, 1.2, 1.2, 0.6, 0.0, 0.1,
+    0.8, 1.1, 0.8, 0.9, 0.0, 0.0, 0.0, 0.7, 0.5,
+    0.4, 0.9, 0.2], dtype=np.float32)
+reemission_cdf = reemission_spec.cumsum()
+reemission_cdf /= reemission_cdf[-1]
+tpb.set('reemission_cdf', reemission_cdf, wavelengths=reemisson_wavelengths)
+
+# Photocathode
+photocathode = Surface('photocathode')
+# Wavelength dependence of R5912-02-MOD PMT, rescaled for higher efficiency
+photocathode.set('detect',
+        np.array([0.0, 0.0, 0.0, 0.0078, 0.039, 0.078, 0.117, 0.1404,
+        0.156, 0.195, 0.2106, 0.195, 0.156, 0.1326, 0.078, 0.039, 0.0195, 0.0078,
+        0.00078, 0.0, 0.0], dtype=np.float32) * 0.25 / 0.16,
+        wavelengths=np.array([60, 200, 260, 270, 280, 285, 290, 300, 310, 330, 370, 420,
+        475, 500, 530, 570, 600, 630, 670, 700, 800], dtype=np.float32))
+
+# Stainless steel: 25% reflectivity below 200 nm, and 50% above.
+# Equal split between specular and diffuse
+steel = Surface("steel")
+steel.set('specular_reflect', [0.125, 0.125, 0.25, 0.25],
+          wavelengths=basic_wavelengths)
+steel.set('diffuse_reflect',  [0.125, 0.125, 0.25, 0.25],
+          wavelengths=basic_wavelengths)
+steel.set('absorb', 1.0 - steel.specular_reflect[:, 1] - steel.diffuse_reflect[:, 1],
+          wavelengths=basic_wavelengths)
+
+# Aluminum (extremely high quality evaporated film)
+aluminum = Surface("aluminum")
+aluminum.set('specular_reflect',
+        np.array([1.0e-3, 1.0e-3, 75.91648e-2, 77.98801e-2, 80.05955e-2,
+        82.28502e-2,  84.36296e-2, 86.2769e-2,  88.66637e-2,  90.58671e-2,
+        92.42781e-2,  94.26708e-2,  95.5589e-2,  96.8507e-2,  97.74627e-2,
+        98.56623e-2,  98.91068e-2,  99.01829e-2,  99.04481e-2,  99.07316e-2,
+        99.09877e-2,  99.05061e-2,  98.91864e-2,  98.86408e-2,  98.65376e-2,
+        98.36511e-2,  97.75578e-2,  97.22937e-2,  96.62096e-2,  96.09272e-2,
+        95.48705e-2,  95.11823e-2,  94.58907e-2,  93.98158e-2,  93.29575e-2,
+        92.60717e-2,  91.76283e-2,  90.91666e-2,  89.75349e-2,  88.66682e-2],
+        dtype=np.float32),
+        wavelengths=np.array([60., 317., 317.56338, 323.98899, 330.41461, 334.25796,
+        344.31949, 351.2497, 358.22436, 368.79049, 379.34921, 388.86909, 402.49245,
+        416.1158, 429.70211, 445.35868, 460.9708, 477.08011, 492.14318, 508.24508,
+        522.78873, 540.44147, 555.48972, 569.50654, 585.06681, 601.13908,
+        615.10404, 631.15408, 645.63846, 660.64968, 676.69231, 692.23775,
+        706.72954, 721.73334, 737.24915, 751.20671, 766.7077, 781.16986,
+        795.60239, 808.48407], dtype=np.float32))

File examples/ly.py

+import sys
+
+abs_length = 1000.0 * int(sys.argv[1])
+from chroma_pds.detector import liquid_Ar, lar5  # , lar5_mask, lar5_mask_double
+
+# Rescale liquid argon absorption length
+liquid_Ar.absorption_length[:2, 1] = abs_length
+print 'liquid argon absorption length (mm):', liquid_Ar.absorption_length
+
+from chroma.loader import create_geometry_from_obj
+from chroma.sim import Simulation
+from chroma.event import *
+from chroma.sample import uniform_sphere
+from chroma.transform import normalize
+import numpy as np
+
+
+def radical_inverse(n, base):
+    val = 0
+    invBase = 1.0 / base
+    invBi = invBase
+
+    while n > 0:
+        d_i = n % base
+        val += d_i * invBi
+        n /= base
+        invBi *= invBase
+    return val
+
+def pick_point_2d(start=0, end=100):
+    for i in xrange(start, end):
+        x, y = np.array([radical_inverse(i, 2), radical_inverse(i, 3)])
+        i += 1
+        yield x, y
+
+
+def pbomb_generator(nphotons, wavelength, pos, dx):
+    # This part doesn't change event to event
+    photon_pos = np.zeros(shape=(nphotons, 3), dtype=np.float32)
+    photon_pos += np.array(pos)
+    photon_pos += np.random.uniform(0.0, 1.0, nphotons)[:, np.newaxis] \
+        * np.asarray(dx)[np.newaxis, :]
+
+    photon_wavelengths = np.empty(shape=nphotons, dtype=np.float32)
+    photon_wavelengths.fill(wavelength)
+
+    photon_t0 = np.zeros_like(photon_wavelengths)
+    photon_dir = uniform_sphere(nphotons).astype(np.float32)
+    photon_pol = normalize(np.cross(uniform_sphere(nphotons), photon_dir)).astype(np.float32)
+    # This part does change
+    while True:
+        yield Photons(pos=photon_pos, dir=photon_dir, pol=photon_pol,
+                      wavelengths=photon_wavelengths, t=photon_t0)
+
+
+def pmt_jitter_scint_dist(fprompt=0.3, fast=6.0, slow=1590.0,
+                          pmt_res=1.5,
+                          xmin=-10, xmax=16000.0, bin_size=0.5):
+    nbins = (xmax-xmin) / bin_size
+    nrounds = 100
+    nphotons = 100000
+
+    total_hist = np.zeros(shape=nbins, dtype=float)
+
+    for i in xrange(nrounds):
+        time_const = np.empty(shape=nphotons, dtype=np.float32)
+        rand_choice = np.random.uniform(0.0, 1.0, nphotons) < fprompt
+        time_const[rand_choice] = fast
+        time_const[~rand_choice] = slow
+        samples = -np.log(np.random.uniform(0.0, 1.0, nphotons)) * time_const
+        samples += np.random.normal(scale=pmt_res, size=samples.shape[0])
+        hist, edges = np.histogram(samples, bins=nbins, range=(xmin, xmax))
+        total_hist += hist
+
+    # All edge arrays will be the same, so we return the last one
+    return edges, total_hist
+
+
+def main(argv):
+    detector = create_geometry_from_obj(lar5())
+    time_edges, time_pdf = pmt_jitter_scint_dist(fprompt=0.3, pmt_res=0.34)
+    detector.set_time_dist(time_edges, time_pdf)
+
+    sim = Simulation(detector, geant4_processes=0)
+    wavelength = 128.0
+    energy = 100
+    nphotons = energy * 19800.0
+
+    with open('ly_scan_newtpb_' + argv[1] + '.txt', 'w') as f:
+        for x, y in pick_point_2d(0, 100000):
+            x *= 5000.0
+            y *= 5000.0
+            pos = (x, y, 0)
+            dx = (0, 0, 7000.0)
+            gen = pbomb_generator(nphotons, wavelength, pos, dx)
+            event = sim.simulate(gen, max_steps=100).next()
+            q = event.channels.q.sum()
+            print x, y, q / energy
+            print >>f, x, y, q / energy
+            f.flush()
+
+if __name__ == '__main__':
+    main(sys.argv)

File examples/plot_ly.py

+import numpy as np
+import ROOT
+import sys
+
+c1 = ROOT.TCanvas('c1', '', 800, 700)
+c1.SetRightMargin(0.15)
+ROOT.gStyle.SetOptStat(0)
+
+scale_factor = 1.0 # 19800.0 / 28000.0
+
+while True:
+    data = np.loadtxt(sys.argv[1])
+    data[:,:2] /= 10.0 # mm -> cm
+    data[:,2] *= scale_factor
+    print 'Rows:', len(data)
+    print 'Average light yield:', data[:,2].mean()
+    nbins = int(sys.argv[2])
+
+
+    h = ROOT.TProfile2D('h', ';Z [beam axis] (cm); X [drift axis] (cm)', 
+                        nbins, data[:,0].min(), data[:,0].max(),
+                        nbins, data[:,1].min(), data[:,1].max())
+    h.SetTitleOffset(1.2, 'Y')
+    for row in data:
+        h.Fill(*row)
+    h.SetMinimum(0.0)
+    h.SetMaximum(1.3)
+    h.Draw('COLZ')
+    c1.Update()
+    c1.Print(sys.argv[3])
+    raw_input()

File examples/steps.py

+from chroma.loader import load_geometry_from_string
+from chroma.sim import Simulation
+from chroma.event import *
+from chroma.sample import uniform_sphere
+from chroma.transform import normalize
+import numpy as np
+import time
+
+
+import sys
+
+
+import pdb
+pdb.set_trace()
+
+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
+    # device, so we call the default hook
+        sys.__excepthook__(type, value, tb)
+    else:
+        import traceback, pdb
+        # we are NOT in interactive mode, print the exception...
+        traceback.print_exception(type, value, tb)
+        print
+        # ...then start the debugger in post-mortem mode.
+        # pdb.pm() # deprecated
+        pdb.post_mortem(tb) # more "modern"
+
+sys.excepthook = info
+
+
+stop_mask = NO_HIT | BULK_ABSORB | SURFACE_ABSORB | SURFACE_DETECT | NAN_ABORT
+
+
+def pbomb_generator(nphotons, wavelength, pos, dx):
+    # This part doesn't change event to event
+    photon_pos = np.zeros(shape=(nphotons, 3), dtype=np.float32)
+    photon_pos += np.array(pos)
+    photon_pos += np.random.uniform(0.0, 1.0, nphotons)[:, np.newaxis] \
+        * np.asarray(dx)[np.newaxis, :]
+
+    photon_wavelengths = np.empty(shape=nphotons, dtype=np.float32)
+    photon_wavelengths.fill(wavelength)
+
+    photon_t0 = np.zeros_like(photon_wavelengths)
+    photon_dir = uniform_sphere(nphotons).astype(np.float32)
+    photon_pol = normalize(np.cross(uniform_sphere(nphotons), photon_dir)).astype(np.float32)
+    # This part does change
+    while True:
+        yield Photons(pos=photon_pos, dir=photon_dir, pol=photon_pol,
+                      wavelengths=photon_wavelengths, t=photon_t0)
+
+mapping = [
+    ['stopped', stop_mask],
+    ['no hit', NO_HIT],
+    ['bulk absorb', BULK_ABSORB],
+    ['surface detect', SURFACE_DETECT],
+    ['surface absorb', SURFACE_ABSORB],
+    ['rayleigh scatter', RAYLEIGH_SCATTER],
+    ['reflect diffuse', REFLECT_DIFFUSE],
+    ['reflect specular', REFLECT_SPECULAR],
+    ['surface reemit', SURFACE_REEMIT],
+    ['surface transmit', SURFACE_TRANSMIT],
+    ['bulk reemit', BULK_REEMIT],
+    ['nan abort', NAN_ABORT]]
+
+
+def print_photon_stats(photons):
+    flags = photons.flags
+    print '%20s:  %d' % ('nphotons', len(flags))
+    print '%20s:  %f' % ('mean lambda',  photons.wavelengths.mean())
+    for label, mask in mapping:
+        selected = (flags & mask) > 0
+        nselected = np.count_nonzero(selected)
+        if nselected > 0:
+            avg_lambda = photons.wavelengths[selected].mean()
+            avg_radius = ((photons.pos[selected] ** 2).sum(axis=1) ** 0.5).mean()
+        else:
+            avg_lambda = 0.0
+            avg_radius = 0.0
+
+        print '%20s:  %d (%1.0f nm, r=%1.1f)' % (label, np.count_nonzero(selected), avg_lambda, avg_radius)
+
+
+geo = load_geometry_from_string('@chroma_pds.detector.lar5')
+
+sim = Simulation(geo, geant4_processes=0)
+
+numpe = []
+
+orig_photons = pbomb_generator(19800 * 10, 128.0,
+                               pos=(-2500, -2500, 0),
+                               dx=(0, 0, -4710.0)).next()
+print_photon_stats(orig_photons)
+photons = orig_photons
+
+for i in xrange(1000):
+    start = time.time()
+
+    event = sim.simulate(photons, keep_photons_beg=True, keep_photons_end=True, max_steps=1).next()
+
+    stop = time.time()
+
+    print_photon_stats(event.photons_end)
+    print
+
+    photons = event.photons_end
+
+    print '(%d) Time: %f' % (i, stop - start)
+import distribute_setup
+distribute_setup.use_setuptools()
+from setuptools import setup, find_packages
+
+setup(
+    name = 'Chroma-PDS',
+    version = '0.1',
+    packages = find_packages(),
+    include_package_data=True,
+
+    install_requires = ['chroma'],
+    test_suite = 'nose.collector',
+)