#!/usr/bin/env python

# Post-processing script that finds events and particle tracks in Geant4
# data. Reads generic Geant4 FITS pixel lists generated by pre-processing
# script 'convet_to_fits_*.py'. Input is a list of raw pixel FITS files,
# including the relative path if needed. The invocation should be of the
# form:
#
#   find_events.py path/rawpix_<runid1>.fits [ path/rawpix_<runid2>.fits ]
#
# where <runid> is a sequence of numbers indicating a unique Geant4 run.
# The script will process each FITS file in turn, sorting primaries into
# radomly selected frames of specified length, and then finding and
# characterizing events and particle tracks ('blobs'). It will write out
# a single FITS file with four extensions:
#
#   path/geant_events_<runid>.fits
#   EVENTS
#       List of events, which are identified as 3x3 islands around
#       a local maximum. 
#   FRAMES
#       Summary information for each frame.
#   BLOBS
#       List of particle tracks, which are contiguous islands of 
#       5 or more pixels (set by npixthresh), or any number of pixels
#       including at least one pixel above the MIP threshold.
#   PIXELS
#       List of all pixels with signal. Basically a copy of the
#       input but with frame and blob info added.
#
# EDM Wed Jul 14 11:41:56 EDT 2021 
# Now handles input files that are gzipped (or at least have .gz
# extension).
#
# EDM Tue Jun  8 16:47:57 EDT 2021
# Removed 'wfits' routine, and now write single output FITS table with 
# four extensions for events, frames, blobs, and pixels.
#
# EDM Wed Mar 31 11:15:02 EDT 2021
# Added primary particle energy PRIMENERGY to output FTIS pixel table.
# Also fixed the fact that PARTYPE, SECPARTYPE, and PRIMTYPE were not
# being populated correctly.
#
# EDM Wed Mar 31 10:38:30 EDT 2021
# Added some elapsed time logging for benchmarking. Also added 'flush=True'
# to most of the print statements so they pop out in the supercloud log
# files as they go.
#
# EDM Tue Mar  9 13:11:35 EST 2021
# Added SECPARTYPE, which records the particle type of the particle that
# entered the detector and ultimately led to the energy deposition. This
# can be the primary or an external secondary, and like PARTYPE it is
# bit-wise added so multiple particle types can be recorded (but not
# multiples of the same particle). PARTYPE encodes final particles, and
# in may cases includes electrons ionized by the proton within the 
# detector and immediately dumping their energy.
#
# EDM Thu Jan 14 14:44:59 EST 2021
# To replicate OU processing, make sure to use the correct values
# of act[x,y]_[min,max], which are commented out below.
#
# EDM Fri Jan  8 13:50:08 EST 2021
# Changed input column EDEP to ENERGY. Probably will break stuff.
#
# EDM Tue Dec 22 08:52:12 EST 2020
# Updated output FITS table data types to reduce file size.  Also
# eliminated EVTTYPE, since it isn't useful and similar filtering can be
# done using other columns. Removed OU-specific hard-wired coordinates in
# event filtering. PTYPE->PARTYPE to eliminate confusion with primary type,
# and added PRIMTYPE to output pixel and event lists. Minor code
# rearrangements.
#
# EDM Mon Dec 14 11:04:27 EST 2020
# First generic version that works. Adapted from 'doit.pl' and 'doit.py'
# scripts that run on the Open University Geant4 data.

import os, sys, glob, re, datetime
import numpy as np
from astropy.table import Table
from astropy.io import fits
from skimage import measure

num_args = len(sys.argv) - 1    
if num_args < 1:
    sys.exit(f"Usage: {sys.argv[0]} <path/to/raw pixel FITS file> ...")

# functions
def match(regex: str, string: str):
    return re.compile(regex).search(string)

def indexND(array, ind, value=None):
    if value is None:
        return array[ind[1],ind[0]]     
    array[ind[1],ind[0]] = value        

def which(condition):
    if condition.ndim != 1:
        condition = condition.flat
    return np.where(condition)[0]

def whichND(condition):
    return np.array(np.where(condition)[::-1])

def which_both(condition):
    return which(condition), which(condition == False)

# define defaults
texp_frame = .005   # frame exposure time in sec
evtth = .1          # 100 eV for WFI Geant4 simulations
splitth = evtth    # same as evtth for WFI Geant4 simulations
npixthresh = 5      # minimum number of pixels in a blob
mipthresh = 15.     # minimum ionizing particle threshold in keV
clip_energy = 22.   # maximum pixel value reported by WFI, in keV
skip_writes = -1    # writes FITS images for every skip_writes primary; 
                    # set to -1 to turn off writes
evperchan = 1000.   # why not? PHA here is really PI
mipthresh_chan = mipthresh * 1000. / evperchan  # minimum ionizing particle threshold in PHA units
spec_maxkev = 100.
numchans = int((spec_maxkev*1000.) / evperchan)
gain_intercept = 0. # use this for Geant4 data
gain_slope = 1000.  # use this for Geant4 data (pixel PH units are keV)
#gain_intercepts = (0, 0, 0, 0)         # in ADU
#gain_slopes = (1000, 1000, 1000, 1000) # in eV/ADU
# rate and frame defaults
proton_flux = 4.1 * 7./5.   # protons/s/cm2; 7/5 accounts for alphas, etc.
detector_area = 4. * (130.e-4 * 512.)**2 # detector area in cm2, 
                                         # assuming 130 um pixels

# conversion from ACIS 256 flt grades to EPIC-pn pattern
# see 'pnpatterns.png' for the patterns
# they include 0 (singles), 1-4 (doubles), 5-8 (triples), 9132 (triples)
epicpn_pattern_table = np.array([
     0, 13,  3, 13, 13, 13, 13, 13,  4, 13,  8, 12, 13, 13, 13, 13,
     2, 13,  7, 13, 13, 13, 11, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
     1, 13, 13, 13, 13, 13, 13, 13,  5, 13, 13, 13, 13, 13, 13, 13,
     6, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13,  9, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    10, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
    13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13 
])

# dictionary of particle types
# these are the same codes as the PDG in MPE data
# and are used for both primary type and particle
# (primary or secondary) that deposits energy;
# can be bit-wise added for multiple secondaries
# depositing in the same pixel for the same primary
ptypes = { 'proton' : 1,
             'e-' : 2,
             'gamma' : 4,
             'e+' : 8,
             'pi-' : 16,
             'pi+' : 32,
             'neutron' : 64,
             'alpha' : 128 }

# initialize rng
rng = np.random.RandomState(1234)

# temporary variables
x, y = 0, 0
x_offset = 0
y_offset = 0

# use these for MIT data
actx_min = 59
actx_max = 1143
acty_min = 52
acty_max = 1121

# use these for OU data for now
#actx_min = 0
#actx_max = 1026
#acty_min = 0
#acty_max = 1026

date_start = datetime.datetime.now().astimezone()
date_string = date_start.strftime("%a %d %b %Y %I:%M:%S %p %Z").rstrip()
print("############################################################", flush=True)
print(f"### Started {sys.argv[0]} on {date_string}", flush=True)
then = date_start

#######################################
# Main loop.
# Step through input raw pixel FITS files.
# For each one, create a frame per primary,
# find blobs and MIPs, and then find events in that frame.

for filename in sys.argv[1:] :

    # parse the filename arg to see if it makes sense
    if not re.search('.*rawpix_[0-9]+\.fits[\.gz]*$', filename) : 
        print(f'Error: file {filename} does not look like a rawpix FITS')
        continue
    # test if the file exists
    if not os.path.isfile(filename) :
        # of course, it could be gzipped
        if not os.path.isfile(f'{filename}.gz') :
            print(f'Error reading file {filename}')
            continue
        else :
            filename = f'{filename}.gz'

    infile = filename
    path = os.path.dirname(filename)
    this_run = int( re.sub(r'.*rawpix_([0-9]+)\.fits.*', r'\1', infile) )
    print(f"### Processing {infile}.", flush=True)

    # read in raw pixel energy deposition FITS file
    # some of this will be written out the pixel list (which will
    # also have frame info)
    rawpix = Table.read(infile)
    pix_actx = rawpix['ACTX']
    pix_acty = rawpix['ACTY']
    pix_energy = rawpix['ENERGY']
    pix_runid = rawpix['RUNID']
    pix_detid = rawpix['DETID']
    pix_primid = rawpix['PRIMID']
    pix_partype = rawpix['PARTYPE']
    pix_secpartype = rawpix['SECPARTYPE']
    pix_primtype = rawpix['PRIMTYPE']
    pix_primenergy = rawpix['PRIMENERGY']
    pix_frame = np.zeros(pix_actx.size, dtype=int)
    pix_blobid = np.zeros(pix_actx.size, dtype=int)

    # get header info
    sphere_radius = rawpix.meta['SPH_RAD']
    numprims_gen = rawpix.meta['NPRI_GEN']
    numprims_interact = rawpix.meta['NPRI_INT']

    # get number of pixels with signal; this is the maximum number
    # we'll use in allocating pixel, event, and blob arrays
    numpix_with_signal = pix_actx.size

    # get unique primids and number of primaries that interact
    uniq_eid = np.unique(pix_primid)
    numprimaries = uniq_eid.size

    # total exposure time in sec for this run
    texp_run = numprims_gen/3.14159/proton_flux/(sphere_radius**2)
    # average num. primaries per frame
    mu_per_frame = numprims_gen * texp_frame / texp_run 

    print(f"### There are {numprims_gen} primaries in this run.", flush=True)
    print(f"### The run exposure time is {texp_run} sec.", flush=True)
    print(f"### The frame exposure time is {texp_frame} sec", flush=True)
    print(f"### for an expected mean of {mu_per_frame} primaries per frame.", flush=True)

    # initialize event piddles, which will be written out or used later
    evt_runid = np.zeros(numpix_with_signal, dtype=np.uint64)
    evt_detid = np.zeros(numpix_with_signal, dtype=np.uint8)
    evt_primid = np.zeros(numpix_with_signal, dtype=np.uint32)
    evt_actx = np.zeros(numpix_with_signal, dtype=np.uint16)
    evt_acty = np.zeros(numpix_with_signal, dtype=np.uint16)
    evt_phas = np.zeros((numpix_with_signal,25), dtype=np.single)
    evt_pha = np.zeros(numpix_with_signal, dtype=np.single)
    evt_partype = np.zeros(numpix_with_signal, dtype=np.uint8)
    evt_secpartype = np.zeros(numpix_with_signal, dtype=np.uint8)
    evt_primtype = np.zeros(numpix_with_signal, dtype=np.uint8)
    evt_primenergy = np.zeros(numpix_with_signal, dtype=np.single)   # in GeV
    evt_energy = np.zeros(numpix_with_signal, dtype = np.single)     # in keV
    evt_blobdist = np.zeros(numpix_with_signal, dtype=np.single)
    evt_mipdist = np.zeros(numpix_with_signal, dtype=np.single)
    evt_pattern = np.zeros(numpix_with_signal, dtype=np.uint8)
    evt_vfaint = np.zeros(numpix_with_signal, dtype=np.uint8)
    # assign frames and times to each primary
    # to start, assume mean of one primary per second
    #evt_time = np.zeros(numpix_with_signal, dtype=float)
    evt_frame = np.zeros(numpix_with_signal, dtype=np.int32)
    pix_time = np.zeros(numpix_with_signal, dtype=float)
    pix_frame = np.zeros(numpix_with_signal, dtype=int)

    # initialize structures to hold the secondary particle columns
    # piddles for the numeric columns; these are enough for now
    #run = np.zeros(0, dtype=int)        # Geant4 run (* in *_detector[0123])
    #detector = np.zeros(0, dtype=int)   # WFI detector (? in *_detector?)
    #eid = np.zeros(0, dtype=int)        # primary ID
    #particleid = np.zeros(0, dtype=int) # interacting particle ID
    #parentid = np.zeros(0, dtype=int)   # don't really need probably

    # initialize piddles to hold blob-specific things to go in FITS table
    blob_frame = np.zeros(numpix_with_signal, dtype=int)
    blob_blobid = np.zeros(numpix_with_signal, dtype=int)
    blob_cenx = np.zeros(numpix_with_signal, dtype=float)
    blob_ceny = np.zeros(numpix_with_signal, dtype=float)
    blob_cenxcl = np.zeros(numpix_with_signal, dtype=float)
    blob_cenycl = np.zeros(numpix_with_signal, dtype=float)
    blob_npix = np.zeros(numpix_with_signal, dtype=int)
    blob_energy = np.zeros(numpix_with_signal, dtype=float)
    blob_energycl = np.zeros(numpix_with_signal, dtype=float)

    # initialize things for the running frames which we will
    # randomly populate
    # frame settings
    # we know there are $numprims_gen, so generate enough
    # random frames to hold them
    framedist = rng.poisson(mu_per_frame, int(2*numprims_gen/mu_per_frame))

    cumframedist = framedist.cumsum()
    # get the total number of frames needed to capture all the primaries; 
    # will write this to FITS header so we can combine runs
    numtotframes = which(cumframedist >= numprims_gen)[0] + 1
    # this is wrong, because it will remove the last bin which we need
    # it's also unnecessary
    #cumframedist = cumframedist[cumframedist <= numprims_gen]

    # initialize piddles to hold frame-specific things to go in FITS table
    frame_frame = np.zeros(numtotframes, dtype=int)
    frame_time = np.zeros(numtotframes, dtype=float)
    frame_runid = np.zeros(numtotframes, dtype=int)
    frame_npix = np.zeros(numtotframes, dtype=int)
    frame_npixmip = np.zeros(numtotframes, dtype=int)
    frame_nevt = np.zeros(numtotframes, dtype=int)
    frame_nevtgood = np.zeros(numtotframes, dtype=int)
    frame_nevt27 = np.zeros(numtotframes, dtype=int)
    frame_nblob = np.zeros(numtotframes, dtype=int)
    frame_nprim = np.zeros(numtotframes, dtype=int)

    # running variables
    numevents = 0
    numtotblobs = 0

    # loop through unique primaries to sort them into frames
    # also figure out the number of primaries with signal in 
    # multiple quadrants, just as a diagnostic
    num_in_different_quadrants = [0, 0, 0, 0]
    """alternate
    for primary in uniq_eid:
    """
    for i in range(numprimaries):
        primary = uniq_eid[i]
        indx = which(pix_primid==primary)
        #print("#####################")
        #print(f"Doing primary {primary}")
        #print(f"indx: {indx}")
        #print(f"pix_primid: {pix_primid[indx]}")
        #print(f"OLD pix_frame {pix_frame[indx]}")
        # assign each primary to a frame
        # first get the frame ID (indexed starts at 0)
        #print(cumframedist[-1])
        frame = which(cumframedist >= primary)[0]
        #print(f"THIS IS FRAME {frame}")
        # then set the primary and pixel piddles
        pix_frame[indx] = frame
        #print(f"NEW pix_frame pix_frame[indx]")
        num_quadrants = np.unique(pix_detid[indx]).size
        num_in_different_quadrants[num_quadrants-1] += 1

    print(f"### Run {this_run}: number of primaries in 1 2 3 4 quads: {num_in_different_quadrants}", flush=True)
    # min and max X,Y values

    # figure out the unique frames that are populated by 
    # primaries that interacted
    frame_frame = np.unique(pix_frame)
    numframes = frame_frame.size
    # now we can allocate frame piddles since we know how many there are
    frame_runid = np.full(numframes, this_run)
    frame_time = np.zeros(numframes, dtype=float)
    frame_npix = np.zeros(numframes)
    frame_npixmip = np.zeros(numframes)
    frame_nevt = np.zeros(numframes)
    frame_nevtgood = np.zeros(numframes)
    frame_nevt27 = np.zeros(numframes)
    frame_nblob = np.zeros(numframes)
    frame_nprim = np.zeros(numframes)

    pct_interact = 100. * numprimaries / numprims_gen
    print(f"### Run {this_run}: generated {numtotframes} total frames,", flush=True)
    print(f"### of which {numframes} frames with the {numprimaries}", flush=True)
    print(f"### interacting primaries will be written.", flush=True)
    print(f"### {numprims_gen} total primaries were simulated.", flush=True)
    print(f"### {numprimaries} or {pct_interact}% of these produced a WFI interaction.", flush=True)

    # loop through frames and make a raw frame for each

    for i in range(numframes):
        # set the frame ID
        this_frame = frame_frame[i]

        # set the frame time
        frame_time[i] = this_frame * texp_frame

        # keep us updated
        if i % 100 == 0:
            now = datetime.datetime.now().astimezone()
            print(f"### Run {this_run}: done {i} of {numframes} frames, {(now-then).total_seconds()} / {(now-date_start).total_seconds()} sec elapsed.", flush=True)
            then = now
            pass

        # are we writing out?
        if skip_writes > 0 and i % skip_writes == 0:
            writeit = 1
        else:
            writeit = 0

        #############################################
        # make a raw frame
        #x = np.zeros(0, float)
        #y = np.zeros(0, float)
        #en = np.zeros(0, float)
        #pixpartype = np.zeros(0, float)
        pixel_indx = which(pix_frame==this_frame)
        x = np.copy(pix_actx[pixel_indx])
        y = np.copy(pix_acty[pixel_indx])
        en = np.copy(pix_energy[pixel_indx])
        pixpartype = np.copy(pix_partype[pixel_indx])
        pixsecpartype = np.copy(pix_secpartype[pixel_indx])
        pixprimtype = np.copy(pix_primtype[pixel_indx])
        pixprimenergy = np.copy(pix_primenergy[pixel_indx])
        pixprimid = np.copy(pix_primid[pixel_indx])
        # we want to populate this so don't sever it
        """comments..."""
        xoff = np.amin(x) + x_offset - 2
        yoff = np.amin(y) + y_offset - 2
        x -= np.amin(x)
        x += 2
        y -= np.amin(y)
        y += 2
        xdim = np.amax(x) + 3
        ydim = np.amax(y) + 3

        #line 542

        # img is the energy (pulseheight) image
        # partypeimg is an image encoding the final particle type(s) responsible
        # secpartypeimg is an image encoding the secondary particle type(s) responsible
        # primidimg is an image encoding the primary responsible
        # for each pixel
        img = np.zeros((ydim, xdim), float) 
        partypeimg = np.zeros((ydim, xdim), float)
        secpartypeimg = np.zeros((ydim, xdim), float)
        primtypeimg = np.zeros((ydim, xdim), float)
        primenergyimg = np.zeros((ydim, xdim), float)
        primidimg = np.zeros((ydim, xdim), float)
        img_xvals = np.fromfunction(lambda x,y: y, (ydim, xdim), dtype = int)
        img_yvals = np.fromfunction(lambda x,y: x, (ydim, xdim), dtype = int)
        #for j in range(en.size):
        #    img[x[j], y[j]] += en[j]
        #    partypeimg[x[j], y[j]] += en[j]
        # better way to do the above mapping to an image without a loop
        # indexND wants a 2xN piddle, where N is the number of pixels
        coos = np.array((x,y))

        indexND(img, coos, en)
        indexND(partypeimg, coos, pixpartype)
        indexND(secpartypeimg, coos, pixsecpartype)
        indexND(primtypeimg, coos, pixprimtype)
        indexND(primenergyimg, coos, pixprimenergy)
        indexND(primidimg, coos, pixprimid)

        #print(pixprimid, x, y, en, sep='\n')
        #print(img)

        # add to some frame piddles
        frame_npix[i] = pixel_indx.size
        frame_npixmip[i] = which(en >= mipthresh).size
        frame_nprim[i] = np.unique(pixprimid).size

        # write out the frame as FITS image
        if writeit:
            fitsfile = os.path.join(path, f'run{this_run}_frame{this_frame}_img.fits')
            print(f"### Writing raw image {fitsfile}.")
            hdu = fits.PrimaryHDU(img)
            hdu.writeto(fitsfile, overwrite=True)

        #############################################
        # find blobs

        # segment image into blobs
        # original IDL code
        # blobRegions=label_region(phimg * (phimg GE evtthresh), /ulong)
        # NLR PDL 'cc8compt' is equivalentish to IDL 'label_regions'
        blobimg = measure.label(img > evtth, connectivity=2)

        # the first blob is 1, not 0 (which indicates no blob)
        # blobadjust decrements blob IDs when we chuck one, so the
        # IDs are continuous from 1
        blobadjust = 0
        for j in range(1, np.amax(blobimg) + 1):
            indx = which(blobimg == j)
            #indx2d = whichND(blobimg == j) dsn
            # set the blob to zeros and skip it if there are too
            # few elements; if it contains a MIP, it's a good blob
            # irregardless
            if (indx.size < npixthresh and np.amax(img.flat[indx]) < mipthresh):
                blobimg.flat[indx] = 0
                blobadjust += 1
                continue
            blobimg.flat[indx] -= blobadjust
            this_blobindx = numtotblobs + j - blobadjust - 1

            # this is the running blobid which we need to add to
            # blob piddles
            blob_blobid[this_blobindx] = numtotblobs + j - blobadjust
            blob_frame[this_blobindx] = this_frame
            blob_npix[this_blobindx] = indx.size
            # calculate unclipped blob centroid and summed energy
            tmp_en = img.flat[indx]
            tmp_toten = np.sum(tmp_en)
            tmp_wtd_x = img_xvals.flat[indx] * tmp_en
            tmp_wtd_y = img_yvals.flat[indx] * tmp_en
            blob_cenx[this_blobindx] = xoff + tmp_wtd_x.sum() / tmp_toten
            blob_ceny[this_blobindx] = yoff + tmp_wtd_y.sum() / tmp_toten
            blob_energy[this_blobindx] = tmp_toten
            # calculate clipped blob centroid and summed energy
            tmp_en = np.clip(tmp_en, None, clip_energy)
            tmp_toten = tmp_en.sum()
            tmp_wtd_x = img_xvals.flat[indx] * tmp_en
            tmp_wtd_y = img_yvals.flat[indx] * tmp_en
            blob_cenxcl[this_blobindx] = xoff + tmp_wtd_x.sum() / tmp_toten
            blob_cenycl[this_blobindx] = yoff + tmp_wtd_y.sum() / tmp_toten
            blob_energycl[this_blobindx] = tmp_toten

        # record number of blobs in this frame
        frame_nblob[i] = np.amax(blobimg)

        # if we found some blobs, change their IDs so they reflect the
        # running total for this run, increase that running blob number 
        blobimg[blobimg > 0] += numtotblobs
        numtotblobs = np.amax(blobimg) if np.amax(blobimg) > 0 else numtotblobs

        # mark each pixel in a blob
        pix_blobid[pixel_indx] = indexND(blobimg, coos)

        # reset the blobimg to ones where there are blobs, zeros
        # otherwise; this way all the blobs have distance 0 automagically
        blobimg[blobimg > 0] = 1

        #############################################
        # find MIPs
        # set up some piddles
        #mip_xy = np.zeros((2, 0), float);

        # segment image into mips
        mipimg = np.copy(img)
        mipimg[mipimg < mipthresh] = 0

        # reset the mipimg to ones where there are mips, zeros
        # otherwise; this way all the mips have distance 0 automagically
        mipimg[mipimg > 0]  = 1

        #############################################
        # find events
        # based on Bev's 'findev.pro' rev3.1, 2011-01-26
        # get indices and number of pixels with PH above event threshold
        # need both 1D and 2D indices, latter to ease local max finding
        indx_thcross = which(img > evtth)
        indx2d_thcross = whichND(img > evtth)
        num_thcross = indx_thcross.size

        #line 676  
        if num_thcross > 0:
            # get piddles containing X,Y coords of threshold crossings
            evtx = indx2d_thcross[0]
            evty = indx2d_thcross[1]
            # ignore threshold crossings on the border, because we don't want to use 
            # these except for local max comparison; also ignore bad pixels and columns
            # as noted with -1 values in the mask image
            # note these indexes are into the $indx_thcross indexes (so they're indexes
            # of indexes)
            # EDM Wed Feb  6 12:56:33 EST 2019 
            # NB: evtx and evty are coords in the small window, not the full frame, so we
            # need to add on the offset.  Also, we want to ignore events on the borders of
            # all segments, not just the outer edge.  And eliminate 2 pixels from all
            # borders.  The last conditional just ensures the pixel is non-negative.  I
            # guess that's important.
            tmp_evtx = evtx + xoff
            tmp_evty = evty + yoff
            indx_notborder = (tmp_evtx>=actx_min+2) & (tmp_evtx<=actx_max-2) & (tmp_evty>=acty_min+2) & (tmp_evty<=acty_max-2) & ( (tmp_evtx<=actx_min+511-2) | (tmp_evtx>=actx_max-511+2) ) & ( (tmp_evty<=acty_min+511-2) | (tmp_evty>=acty_max-511+2) )
            #indx_border, indx_notborder = which_both(
            #    (tmp_evtx<=1) | ((tmp_evtx>=510) & (tmp_evtx<=516)) | (tmp_evtx>=1025) | 
            #    (tmp_evty<=1) | ((tmp_evty>=510) & (tmp_evty<=516)) | (tmp_evty>=1025) | 
            #    (img[evty,evtx] < 0)
            #)
            num_notborder = indx_notborder.size

            if num_notborder > 0:
                evtx = evtx[indx_notborder]
                evty = evty[indx_notborder]
                indx_thcross = indx_thcross[indx_notborder]
                # find local maxima
                # make a copy of the threshold crossing piddle, which we will
                # subsitute with the largest neighbor value and then compare to the original
                localmax = np.copy(img)
                """comments..."""

                for xi,yi in [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]:             
                    indx = which(img[evty + yi,evtx + xi] > localmax[evty,evtx])
                    localmax.flat[indx_thcross[indx]] = img[evty + yi, evtx + xi].flat[indx]

                # finally compare the original central pixel pulseheight
                # with the maximum neighbor; if the former is greater than
                # or equal to the latter, it's a local maximum
                indx_localmax = which(img.flat[indx_thcross] >= localmax.flat[indx_thcross])
                num_localmax = indx_localmax.size

                #line 748
                if (num_localmax > 0):
                    evtx = evtx[indx_localmax]
                    evty = evty[indx_localmax]
                    indx_thcross = indx_thcross[indx_localmax]
                    # get 3x3 PHAS piddle in correct
                    # 1D order (central pixel first)
                    this_phas = np.zeros((num_localmax, 25), float)

                    phas_order = iter(range(1, 25))
                    this_phas[:,0] = img[evty, evtx]  
                    for yi in range(-1,2):
                        for xi in range(-1,2):
                            if (xi, yi) != (0,0):
                                this_phas[:,next(phas_order)] = img[evty + yi,evtx + xi]              

                    # get outer 5x5 PHAS piddle
                    # first pad the image with zeros all around
                    tmp_img = np.pad(img, 1, mode='constant')
                    # now evtx and evty are too low by one, so correct for that in index
                    for yi in range(-1,4):
                        for xi in range(-1,4):
                            if not(-1 < xi < 3 and -1 < yi < 3):
                                this_phas[:, next(phas_order)] = tmp_img[evty + yi,evtx + xi]

                    # set the particle types for all these events
                    # based on whatever partype produced the local max
                    this_partype = np.copy(partypeimg[evty,evtx])
                    this_secpartype = np.copy(secpartypeimg[evty,evtx])
                    this_primtype = np.copy(primtypeimg[evty,evtx])
                    this_primenergy = np.copy(primenergyimg[evty,evtx])
                    this_primid = np.copy(primidimg[evty,evtx])
                else:
                    continue
            else:
                continue
        else:
            continue

        ## get X,Y coords of any pixels in a blob (don't matter which blob)
        indx_blobs = whichND(blobimg > 0).T

        ## get X,Y coords of any pixels in a mip (don't matter which mip)
        indx_mips = whichND(mipimg > 0).T

        numevents_thisframe = num_localmax
        #print(f"### Found {numevents_thisframe} events, processing.")

        # process the detected events
        # append events to running piddles
        indx = np.arange(numevents,numevents+numevents_thisframe)
        evt_runid[indx] = this_run
        #evt_detid[indx] = 0
        evt_primid[indx] = this_primid
        evt_frame[indx] = this_frame
        evt_actx[indx] = evtx
        evt_acty[indx] = evty
        evt_partype[indx] = this_partype
        evt_secpartype[indx] = this_secpartype
        evt_primtype[indx] = this_primtype
        evt_primenergy[indx] = this_primenergy
        evt_phas[indx,:] = this_phas
        #evt_pha[indx] = np.zeros(numevents_thisframe)
        #evt_partype[indx] = np.array(this_partype, np.uint8)
        #evt_energy[indx] = np.zeros(numevents_thisframe)
        #evt_blobdist[indx] = np.zeros(numevents_thisframe)
        #evt_mipdist[indx] = np.zeros(numevents_thisframe)
        #evt_pattern[indx] = np.zeros(numevents_thisframe, np.uint8)
        #evt_vfaint[indx] = np.zeros(numevents_thisframe, np.uint8)

        # step through all events to determine
        # EPIC-pn pattern, summed PHA, VFAINT flag, minimum distance to a
        # blob and mip.
        for j in range(numevents, numevents+numevents_thisframe):
            # below is already done in event finding and pasted above
            # get X,Y of center pixel
            x = evt_actx[j]
            y = evt_acty[j]
            # below is deprecated, since we've assumed this in event finding
            ## eliminate events on edges so we can use 5x5
            #next if ($x<2 or $x>1021 or $y<2 or $y>1021);

            # get ACIS flt grade; "which" returns the indices of
            # non-central pixels greater than or equal to event threshold,
            # these are used as bits to raise 2 to the power, and summed
            # (this can probably be removed from the loop if I'm smart)
            indx = which(evt_phas[j,1:9] >= splitth)
            fltgrade = pow(2,indx).sum()
            # convert to EPIC-pn pattern (from look-up table)
            evt_pattern[j] = epicpn_pattern_table[fltgrade]

            # sum 3x3 pixels over split threshold to get PHA
            # this is an ACIS PHA, not Suzaku
            # (this can probably be removed from the loop if I'm smart)
            evt_pha[j] = evt_phas[j][evt_phas[j] >= splitth].sum()
            # apply gain correction for this node
            # get the gain parameters from the node (0-3)
            #print(f"{evt_phas[j]} {gain_intercept} {gain_slope} {evt_pha[j]}")
            evt_pha[j] = (evt_pha[j] - gain_intercept) * gain_slope / evperchan

            # convert pha to energy
            # (this can probably be removed from the loop if I'm smart)
            evt_energy[j] = evt_pha[j] * evperchan / 1000
            #print(f"{(evt_phas[j])} {gain_intercept} {gain_slope} {evt_pha[j]} {evt_energy[j]}")

            # perform simple VFAINT filtering; also update the EPIC-pn pattern
            # of doubles, triples, and quads based on it
            # get outer 16 pixels and if any are above split flag it

            noutermost = which(evt_phas[j, 9:25] > splitth).size #maybe don't need which
            if noutermost > 0:
                evt_vfaint[j] = 1
                # EDM Fri Jan  4 11:24:20 EST 2019 
                # Reinstated the 5x5 filtering on PATTERN.
                # EDM Thu May 31 13:46:28 EDT 2018 
                # change to remove 5x5 criterion on EPIC-pn patterns
                # for doubles, triples, quads
                if evt_pattern[j] > 0: 
                    evt_pattern[j] = 13

            # get minimum distance to a blob
            # first find delta X,Y from this event to list of blob pixels
            delta_blobs = indx_blobs - np.array([x,y])
            #print(delta_blobs)
            # square that, sum it, square root it to get the distance
            if delta_blobs.shape[0] > 0:
                evt_blobdist[j] = np.amin(np.sqrt(np.square(delta_blobs).sum(axis = 1)))
            # unless there aren't any blobs, in which case set blobdist to -1
            else:
                evt_blobdist[j] = -1
            #print(f"{evt_blobdist[j]}")

            # get minimum distance to a mip
            # first find delta X,Y from this event to list of mip pixels
            delta_mips = indx_mips - np.array([x,y])
            #print $delta_mips;
            # square that, sum it, square root it to get the distance
            if delta_mips.shape[0] > 0:

                evt_mipdist[j] = np.amin(np.sqrt(np.square(delta_mips).sum(axis = 1)))
            # unless there aren't any mips, in which case set mipdist to -1
            else:
                evt_mipdist[j] = -1

            #print(f"{evt_mipdist[j]}")

            # we really want ACTX,Y in real WFI coords in the event list, 
            # so fix that here; $x and $y are children of $actx and $acty
            # for this event, which is why this works
            evt_actx[j] += xoff
            evt_acty[j] += yoff

            # add info to frame piddles
            frame_nevt[i] += 1
            if evt_energy[j] < mipthresh:
                frame_nevtgood[i] += 1
            if evt_energy[j] >= 2 and evt_energy[j] < 7:
                frame_nevt27[i] += 1


        # increment the number of events and hit the next frame
        numevents += numevents_thisframe
    # done loop through frames

    # set header/metadata for output FITS files
    hdr = fits.Header()
    hdr.set('SPH_RAD', sphere_radius, 'Radius of Geant4 source sphere in cm')
    hdr.set('NPRI_GEN', numprims_gen, 'Number of primaries generated')
    hdr.set('NPRI_INT', numprims_interact, 'Number of primaries producing signal')
    hdr.set('NFRAMES', numtotframes, 'Number of frames in this file')
    hdr.set('NBLOBS', numtotblobs, 'Number of blobs in this file')
    hdr.set('TEXPFRAM', texp_frame, 'Frame exposure time in sec')

    # create an empty primary HDU
    empty_primary = fits.PrimaryHDU(header=hdr)

    # create events HDU, first chopping off the unused parts of the arrays
    col1  = fits.Column(name='FRAME', format='J', array=evt_frame[:numevents])
    col2  = fits.Column(name='DETID', format='B', array=evt_detid[:numevents])
    col3  = fits.Column(name='ACTX', format='I', array=evt_actx[:numevents])
    col4  = fits.Column(name='ACTY', format='I', array=evt_acty[:numevents])
    col5  = fits.Column(name='ENERGY', format='E', array=evt_energy[:numevents])
    col6  = fits.Column(name='BLOBDIST', format='E', array=evt_blobdist[:numevents])
    col7  = fits.Column(name='MIPDIST', format='E', array=evt_mipdist[:numevents])
    col8  = fits.Column(name='PATTERN', format='B', array=evt_pattern[:numevents])
    col9  = fits.Column(name='PARTYPE', format='B', array=evt_partype[:numevents])
    col10 = fits.Column(name='SECPARTYPE', format='B', array=evt_secpartype[:numevents])
    col11 = fits.Column(name='PRIMTYPE', format='B', array=evt_primtype[:numevents])
    col12 = fits.Column(name='PRIMENERGY', format='E', array=evt_primenergy[:numevents])
    col13 = fits.Column(name='VFAINT', format='B', array=evt_vfaint[:numevents])
    col14 = fits.Column(name='PHA', format='E', array=evt_pha[:numevents])
    col15 = fits.Column(name='PHAS', format='25E', array=evt_phas[:numevents])
    col16 = fits.Column(name='RUNID', format='K', array=evt_runid[:numevents])
    col17 = fits.Column(name='PRIMID', format='J', array=evt_primid[:numevents])
    events_hdu = fits.BinTableHDU.from_columns([col1, col2, col3, col4, col5, col6, col7, col8,
        col9, col10, col11, col12, col13, col14, col15, col16, col17], name='EVENTS', header=hdr)

    # create pixel HDU, after sorting by frame
    pix_actx += x_offset
    pix_acty += y_offset
    indx_sorted = np.argsort(pix_frame)
    col1 = fits.Column(name='FRAME', format='J', array=pix_frame[indx_sorted])
    col2 = fits.Column(name='DETID', format='B', array=pix_detid[indx_sorted])
    col3 = fits.Column(name='ACTX', format='I', array=pix_actx[indx_sorted])
    col4 = fits.Column(name='ACTY', format='I', array=pix_acty[indx_sorted])
    col5 = fits.Column(name='ENERGY', format='E', array=pix_energy[indx_sorted])
    col6 = fits.Column(name='BLOBID', format='K', array=pix_blobid[indx_sorted])
    col7 = fits.Column(name='PARTYPE', format='B', array=pix_partype[indx_sorted])
    col8 = fits.Column(name='SECPARTYPE', format='B', array=pix_secpartype[indx_sorted])
    col9 = fits.Column(name='PRIMTYPE', format='B', array=pix_primtype[indx_sorted])
    col10 = fits.Column(name='PRIMENERGY', format='E', array=pix_primenergy[indx_sorted]) 
    col11 = fits.Column(name='RUNID', format='K', array=pix_runid[indx_sorted])      
    col12 = fits.Column(name='PRIMID', format='J', array=pix_primid[indx_sorted])     
    pix_hdu = fits.BinTableHDU.from_columns([col1, col2, col3, col4, col5, col6, col7, col8,
        col9, col10, col11, col12], name='PIXELS', header=hdr)

    # create frames HDU, after sorting by frame
    indx_sorted = np.argsort(frame_frame)
    col1 = fits.Column(name='FRAME', format='J', array=frame_frame[indx_sorted])
    col2 = fits.Column(name='TIME', format='E', array=frame_time[indx_sorted])
    col3 = fits.Column(name='NPRIM', format='I', array=frame_nprim[indx_sorted])
    col4 = fits.Column(name='NBLOB', format='I', array=frame_nblob[indx_sorted])
    col5 = fits.Column(name='NEVT', format='I', array=frame_nevt[indx_sorted])
    col6 = fits.Column(name='NEVT27', format='I', array=frame_nevt27[indx_sorted])
    col7 = fits.Column(name='NEVTGOOD', format='I', array=frame_nevtgood[indx_sorted])
    col8 = fits.Column(name='NPIX', format='J', array=frame_npix[indx_sorted])
    col9 = fits.Column(name='NPIXMIP', format='J', array=frame_npixmip[indx_sorted])
    col10 = fits.Column(name='RUNID', format='K', array=frame_runid[indx_sorted])
    frames_hdu = fits.BinTableHDU.from_columns([col1, col2, col3, col4, col5, col6, col7, col8,
        col9, col10], name='FRAMES', header=hdr)

    # create blobs HDU, after chop off the unused parts of the arrays and sorting by blobid
    blob_frame = blob_frame[:numtotblobs]
    blob_blobid = blob_blobid[:numtotblobs]
    blob_cenx = blob_cenx[:numtotblobs]
    blob_ceny = blob_ceny[:numtotblobs]
    blob_cenxcl = blob_cenxcl[:numtotblobs]
    blob_cenycl = blob_cenycl[:numtotblobs]
    blob_npix = blob_npix[:numtotblobs]
    blob_energy = blob_energy[:numtotblobs]
    blob_energycl = blob_energycl[:numtotblobs]
    indx_sorted = np.argsort(blob_blobid)
    col1 = fits.Column(name='FRAME', format='J', array=blob_frame[indx_sorted])
    col2 = fits.Column(name='BLOBID', format='K', array=blob_blobid[indx_sorted])
    col3 = fits.Column(name='CENX', format='E', array=blob_cenx[indx_sorted])
    col4 = fits.Column(name='CENXCL', format='E', array=blob_cenxcl[indx_sorted])
    col5 = fits.Column(name='CENY', format='E', array=blob_ceny[indx_sorted])
    col6 = fits.Column(name='CENYCL', format='E', array=blob_cenycl[indx_sorted])
    col7 = fits.Column(name='ENERGY', format='E', array=blob_energy[indx_sorted])
    col8 = fits.Column(name='ENERGYCL', format='E', array=blob_energycl[indx_sorted])
    col9 = fits.Column(name='NPIX', format='J', array=blob_npix[indx_sorted])
    blobs_hdu = fits.BinTableHDU.from_columns([col1, col2, col3, col4, col5, col6, col7, col8,
        col9], name='BLOBS', header=hdr)

    # write the FITS table
    print(f"### Writing output FITS file for run {this_run}.")
    outfile = os.path.join(path, f'geant_events_{this_run}.fits')
    hdul = fits.HDUList([empty_primary, events_hdu, frames_hdu, blobs_hdu, pix_hdu])
    hdul.writeto(outfile, output_verify='silentfix', overwrite=True, checksum=True)

    print(f"### Finished run {this_run}.")

# done loop through runs

date_stop = datetime.datetime.now().astimezone()
elapsed = (date_stop-date_start).total_seconds()
date_string = date_stop.strftime("%a %d %b %Y %I:%M:%S %p %Z").rstrip()
print(f"### Finished {sys.argv[0]} on {date_string}, took {elapsed} sec elapsed.", flush=True)
print("############################################################", flush=True)
