def tt_bblock (tstarts, tstops, times, p0=0.05, intersect_with_bins=False):
    """Bayesian Blocks for time-tagged events. Arguments:

    *tstarts*
      Array of input bin start times.
    *tstops*
      Array of input bin stop times.
    *times*
      Array of event arrival times.
    *p0* = 0.05
      Probability of preferring solutions with additional bins.
    *intersect_with_bins* = False
      If true, intersect bblock bins with input bins; can result in more bins
      than bblocks wants; they will have the same rate values.

    Returns a Holder with:

    *counts*
      Number of events in each output block.
    *finalp0*
      Final value of p0, after iteration to minimize `nblocks`.
    *ledges*
      Times of left edges of output blocks.
    *midpoints*
      Times of midpoints of output blocks.
    *nblocks*
      Number of output blocks.
    *ncells*
      Number of input cells/bins.
    *origp0*
      Original value of p0.
    *rates*
      Event rate associated with each block.
    *redges*
      Times of right edges of output blocks.
    *widths*
      Width of each output block.

    Bin start/stop times are best derived from a 1D Voronoi tesselation of the
    event arrival times, with some kind of global observation start/stop time
    setting the extreme edges. Or they can be set from "good time intervals"
    if observations were toggled on or off as in an X-ray telescope.

    If *intersect_with_bins* is True, the true Bayesian Blocks bins (BBBs) are
    intersected with the "good time intervals" (GTIs) defined by the *tstarts*
    and *tstops* variables. One GTI may contain multiple BBBs if the event
    rate appears to change within the GTI, and one BBB may span multiple GTIs
    if the event date does *not* appear to change between the GTIs. The
    intersection will ensure that no BBB intervals cross the edge of a GTI. If
    this would happen, the BBB is split into multiple, partially redundant
    records. Each of these records will have the **same** value for the
    *counts*, *rates*, and *widths* values. However, the *ledges*, *redges*,
    and *midpoints* values will be recalculated. Note that in this mode, it is
    not necessarily true that ``widths = ledges - redges`` as is usually the
    case. When this flag is true, keep in mind that multiple bins are
    therefore *not* necessarily independent statistical samples.

    """
    tstarts = np.asarray (tstarts)
    tstops = np.asarray (tstops)
    times = np.asarray (times)

    if tstarts.size != tstops.size:
        raise ValueError ('must have same number of starts and stops')

    ngti = tstarts.size

    if ngti < 1:
        raise ValueError ('must have at least one goodtime interval')
    if np.any ((tstarts[1:] - tstarts[:-1]) <= 0):
        raise ValueError ('tstarts must be ordered and distinct')
    if np.any ((tstops[1:] - tstops[:-1]) <= 0):
        raise ValueError ('tstops must be ordered and distinct')
    if np.any (tstarts >= tstops):
        raise ValueError ('tstarts must come before tstops')
    if np.any ((times[1:] - times[:-1]) < 0):
        raise ValueError ('times must be ordered')
    if times.min () < tstarts[0]:
        raise ValueError ('no times may be smaller than first tstart')
    if times.max () > tstops[-1]:
        raise ValueError ('no times may be larger than last tstop')
    for i in range (1, ngti):
        if np.where ((times > tstops[i-1]) & (times < tstarts[i]))[0].size:
            raise ValueError ('no times may fall in goodtime gap #%d' % i)
    if p0 < 0 or p0 >= 1.:
        raise ValueError ('p0 must lie within [0, 1)')

    utimes, uidxs = np.unique (times, return_index=True)
    nunique = utimes.size

    counts = np.empty (nunique)
    counts[:-1] = uidxs[1:] - uidxs[:-1]
    counts[-1] = times.size - uidxs[-1]
    assert counts.sum () == times.size

    # we grow these arrays with concats, which will perform badly with lots of
    # GTIs. Not expected to be a big deal.
    widths = np.empty (0)
    ledges = np.empty (0)
    redges = np.empty (0)

    for i in range (ngti):
        tstart, tstop = tstarts[i], tstops[i]

        w = np.where ((utimes >= tstart) & (utimes <= tstop))[0]

        if not w.size:
            # No events during this goodtime! We have to insert a zero-count
            # event block. This may break assumptions within bin_bblock()?

            # j = idx of first event after this GTI
            wafter = np.where (utimes > tstop)[0]
            if wafter.size:
                j = wafter[0]
            else:
                j = utimes.size
            assert j == 0 or np.where (utimes < tstart)[0][-1] == j - 1

            counts = np.concatenate ((counts[:j], [0], counts[j:]))
            widths = np.concatenate ((widths, [tstop - tstart]))
            ledges = np.concatenate ((ledges, [tstart]))
            redges = np.concatenate ((redges, [tstop]))
        else:
            gtutimes = utimes[w]
            midpoints = 0.5 * (gtutimes[1:] + gtutimes[:-1]) # size: n - 1
            gtedges = np.concatenate (([tstart], midpoints, [tstop])) # size: n + 1
            gtwidths = gtedges[1:] - gtedges[:-1] # size: n
            assert gtwidths.sum () == tstop - tstart
            widths = np.concatenate ((widths, gtwidths))
            ledges = np.concatenate ((ledges, gtedges[:-1]))
            redges = np.concatenate ((redges, gtedges[1:]))

    assert counts.size == widths.size
    info = bin_bblock (widths, counts, p0=p0)
    info.ledges = ledges[info.blockstarts]
    # The right edge of the i'th block is the right edge of its rightmost
    # bin, which is the bin before the leftmost bin of the (i+1)'th block:
    info.redges = np.concatenate ((redges[info.blockstarts[1:] - 1], [redges[-1]]))
    info.midpoints = 0.5 * (info.ledges + info.redges)
    del info.blockstarts

    if intersect_with_bins:
        # OK, we now need to intersect the bblock bins with the input bins.
        # This can fracture one bblock bin into multiple ones but shouldn't
        # make any of them disappear, since they're definitionally constrained
        # to contain events.
        #
        # First: sorted list of all timestamps at which *something* changes:
        # either a bblock edge, or a input bin edge. We drop the last entry,
        # giving is a list of left edges of bins in which everything is the
        # same.

        all_times = set(tstarts)
        all_times.update(tstops)
        all_times.update(info.ledges)
        all_times.update(info.redges)
        all_times = np.array(sorted(all_times))[:-1]

        # Now, construct a lookup table of which bblock number each of these
        # bins corresponds to. More than one bin may have the same bblock
        # number, if a GTI change slices a single bblock into more than one
        # piece. We do this in a somewhat non-obvious way since we know that
        # the bblocks completely cover the overall GTI span in order.

        bblock_ids = np.zeros(all_times.size)

        for i in range(1, info.nblocks):
            bblock_ids[all_times >= info.ledges[i]] = i

        # Now, a lookup table of which bins are within a good GTI span. Again,
        # we know that all bins are either entirely in a good GTI or entirely
        # outside, so the logic is simplified but not necessarily obvious.

        good_timeslot = np.zeros(all_times.size, dtype=np.bool)

        for t0, t1 in zip(tstarts, tstops):
            ok = (all_times >= t0) & (all_times < t1)
            good_timeslot[ok] = True

        # Finally, look for contiguous spans that are in a good timeslot *and*
        # have the same underlying bblock number. These are our intersected
        # blocks.

        old_bblock_ids = []
        ledges = []
        redges = []
        cur_bblock_id = -1

        for i in range(all_times.size):
            if bblock_ids[i] != cur_bblock_id or not good_timeslot[i]:
                if cur_bblock_id >= 0:
                    # Ending a previous span.
                    redges.append(all_times[i])
                    cur_bblock_id = -1

                if good_timeslot[i]:
                    # Starting a new span.
                    ledges.append(all_times[i])
                    old_bblock_ids.append(bblock_ids[i])
                    cur_bblock_id = bblock_ids[i]

        if cur_bblock_id >= 0:
            # End the last span.
            redges.append(tstops[-1])

        # Finally, rewrite all of the data as planned.

        old_bblock_ids = np.array(old_bblock_ids, dtype=np.int)
        info.counts = info.counts[old_bblock_ids]
        info.rates = info.rates[old_bblock_ids]
        info.widths = info.widths[old_bblock_ids]

        info.ledges = np.array(ledges)
        info.redges = np.array(redges)
        info.midpoints = 0.5 * (info.ledges + info.redges)
        info.nblocks = info.ledges.size

    return info