| 1 | #!/usr/bin/python -tt
|
|---|
| 2 | #
|
|---|
| 3 | #
|
|---|
| 4 |
|
|---|
| 5 | from pyfact import RawData
|
|---|
| 6 | from drs_spikes import DRSSpikes
|
|---|
| 7 | from fir_filter import CFD
|
|---|
| 8 | from fir_filter import SlidingAverage
|
|---|
| 9 |
|
|---|
| 10 | import os.path
|
|---|
| 11 | import sys
|
|---|
| 12 | import numpy as np
|
|---|
| 13 | try:
|
|---|
| 14 | import h5py
|
|---|
| 15 | except:
|
|---|
| 16 | print 'module h5py needed for this process ... aborting'
|
|---|
| 17 | sys.exit(1)
|
|---|
| 18 |
|
|---|
| 19 | data_filename = 'data/20111017_010.fits.gz'
|
|---|
| 20 | calib_filename = 'data/20111017_006.drs.fits.gz'
|
|---|
| 21 |
|
|---|
| 22 | run = RawData(data_filename, calib_filename, return_dict = True, do_calibration=True)
|
|---|
| 23 | despike = DRSSpikes()
|
|---|
| 24 | smooth = SlidingAverage(7)
|
|---|
| 25 | cfd = CFD()
|
|---|
| 26 |
|
|---|
| 27 | thr = 3
|
|---|
| 28 | filter_delay = 3
|
|---|
| 29 | search_window_size = 12
|
|---|
| 30 | # shortcut
|
|---|
| 31 | sws = search_window_size
|
|---|
| 32 |
|
|---|
| 33 |
|
|---|
| 34 | # we try to determine the:
|
|---|
| 35 | # * Position
|
|---|
| 36 | # * and Height (in two versions: filtered, and unfiltered data)
|
|---|
| 37 | # of dark count peaks in:
|
|---|
| 38 | # a) all events
|
|---|
| 39 | # b) all pixel
|
|---|
| 40 | #
|
|---|
| 41 | # we limit ourselfs to 10 peaks per pixel max!
|
|---|
| 42 | # and 1000 events max!
|
|---|
| 43 | #
|
|---|
| 44 | # So in order to store this stuff, we prepare an n-dim-array
|
|---|
| 45 |
|
|---|
| 46 | PEAKS_MAX = 30
|
|---|
| 47 |
|
|---|
| 48 | # In order not to allocate so much mem, we use hdf5 files for storing this information
|
|---|
| 49 | f = h5py.File('20111017_010-006.h5', 'w')
|
|---|
| 50 |
|
|---|
| 51 | # we write some meta data into the file, so we know afterwards, where this stuff
|
|---|
| 52 | # came from
|
|---|
| 53 | f.attrs['data_file'] = os.path.realpath(data_filename)
|
|---|
| 54 | f.attrs['calib_file'] = os.path.realpath(calib_filename)
|
|---|
| 55 | # this script
|
|---|
| 56 | f.attrs['source_python_script'] = open(sys.argv[0], 'r').read()
|
|---|
| 57 | # SVN revision, should go here ... but of what? this file??
|
|---|
| 58 | f.attrs['settings'] = [
|
|---|
| 59 | 'threshold = 3',
|
|---|
| 60 | 'SlidingAverage = 7 --> Filter delay = 3'
|
|---|
| 61 | 'Search Window Size = 12'
|
|---|
| 62 | 'Maximum Number of Peaks, in data sets = 30'
|
|---|
| 63 | 'Default value of peak_positions = -1'
|
|---|
| 64 | 'Default value of peak_height_unfiltered = -np.inf'
|
|---|
| 65 | 'Default value of peak_height_smoothed = -np.inf'
|
|---|
| 66 | ]
|
|---|
| 67 |
|
|---|
| 68 | result_peak_positions = f.create_dataset('peak_positions',
|
|---|
| 69 | (run.nevents, run.npix, PEAKS_MAX),
|
|---|
| 70 | dtype=np.int16,
|
|---|
| 71 | fillvalue=-1,
|
|---|
| 72 | compression='gzip')
|
|---|
| 73 | result_peak_unfiltered_height = f.create_dataset('peak_height_unfiltered',
|
|---|
| 74 | (run.nevents, run.npix, PEAKS_MAX),
|
|---|
| 75 | dtype=np.float32,
|
|---|
| 76 | fillvalue=-np.inf,
|
|---|
| 77 | compression='gzip')
|
|---|
| 78 | result_peak_smoothed_height = f.create_dataset('peak_height_smoothed',
|
|---|
| 79 | (run.nevents, run.npix, PEAKS_MAX),
|
|---|
| 80 | dtype=np.float32,
|
|---|
| 81 | fillvalue=-np.inf,
|
|---|
| 82 | compression='gzip')
|
|---|
| 83 |
|
|---|
| 84 | try:
|
|---|
| 85 |
|
|---|
| 86 | for event in run:
|
|---|
| 87 | event_id = event['event_id'].value - 1
|
|---|
| 88 |
|
|---|
| 89 | data = event['acal_data']
|
|---|
| 90 | data = despike(data)
|
|---|
| 91 | data_orig = data.copy()
|
|---|
| 92 | data = smooth(data)
|
|---|
| 93 | filtered = cfd(data)
|
|---|
| 94 | filtered = smooth(filtered)
|
|---|
| 95 |
|
|---|
| 96 | # this is a loop over all pixel of this event
|
|---|
| 97 | pixel_id = -1
|
|---|
| 98 | for dat, fil, orig in zip(data, filtered, data_orig):
|
|---|
| 99 | pixel_id += 1
|
|---|
| 100 | print event_id, pixel_id
|
|---|
| 101 |
|
|---|
| 102 | prod = fil[:-1] * fil[1:]
|
|---|
| 103 | cand = np.where( prod <= 0)[0]
|
|---|
| 104 | if len(cand) == 0:
|
|---|
| 105 | continue
|
|---|
| 106 | # zero crossing with rising edge
|
|---|
| 107 | cross = cand[np.where(fil[cand] < 0)[0]]
|
|---|
| 108 | if len(cross) == 0:
|
|---|
| 109 | continue
|
|---|
| 110 | over_thr = cross[np.where(dat[cross-4] > thr)[0]]
|
|---|
| 111 |
|
|---|
| 112 | # Now since we have these values, we will throw away all those,
|
|---|
| 113 | # which are probably on a falling edge of its predecessor
|
|---|
| 114 | dover = np.diff(over_thr)
|
|---|
| 115 | if len(dover) == 0:
|
|---|
| 116 | good = over_thr
|
|---|
| 117 | else:
|
|---|
| 118 | good = []
|
|---|
| 119 | good.append(over_thr[0])
|
|---|
| 120 | for i in range(len(dover)):
|
|---|
| 121 | if dover[-i-1] > 100:
|
|---|
| 122 | good.append(over_thr[-i-1])
|
|---|
| 123 |
|
|---|
| 124 | good = np.array(sorted(good))
|
|---|
| 125 | # these positions, we just found, do not exactly point to the maximum
|
|---|
| 126 | # of a peak, but the peak will be on the left side of it.
|
|---|
| 127 | # we use the smoothed data to find the position of the local maximum
|
|---|
| 128 | # and then stopre this position and the value of both
|
|---|
| 129 | # the smoothed data and the original data.
|
|---|
| 130 |
|
|---|
| 131 | max_pos = np.zeros( good.shape )
|
|---|
| 132 | max_smoothed = np.zeros( good.shape )
|
|---|
| 133 | max_orig = np.zeros( good.shape )
|
|---|
| 134 |
|
|---|
| 135 | for i in range(len(good)):
|
|---|
| 136 | # We search for a local maximum in a window of size 12
|
|---|
| 137 | if len(dat[good[i]-sws:good[i]]) > 0:
|
|---|
| 138 |
|
|---|
| 139 | max_pos[i] = good[i]-sws + np.argmax(dat[good[i]-sws:good[i]])
|
|---|
| 140 | max_smoothed[i] = dat[max_pos[i]]
|
|---|
| 141 | max_orig[i] = orig[max_pos[i]-filter_delay]
|
|---|
| 142 |
|
|---|
| 143 | result_peak_positions[event_id,pixel_id, :len(max_pos)] = max_pos
|
|---|
| 144 | result_peak_unfiltered_height[event_id,pixel_id, :len(max_pos)] =max_orig
|
|---|
| 145 | result_peak_smoothed_height[event_id,pixel_id, :len(max_pos)] = max_smoothed
|
|---|
| 146 |
|
|---|
| 147 |
|
|---|
| 148 | f.close()
|
|---|
| 149 |
|
|---|
| 150 | except:
|
|---|
| 151 | print 'An Exception occured .. closing file'
|
|---|
| 152 | f.close()
|
|---|
| 153 | raise
|
|---|
| 154 |
|
|---|
| 155 |
|
|---|
| 156 |
|
|---|
