| 1 | #!/usr/bin/python -tti
|
|---|
| 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 | from plotters import Plotter
|
|---|
| 11 |
|
|---|
| 12 | import sys
|
|---|
| 13 | import cPickle
|
|---|
| 14 | import matplotlib.pyplot as plt
|
|---|
| 15 | import numpy as np
|
|---|
| 16 |
|
|---|
| 17 | data_filename = 'data/20111017_009.fits.gz'
|
|---|
| 18 | calib_filename = 'data/20111017_006.drs.fits.gz'
|
|---|
| 19 | out_filename = 'test.pkl'
|
|---|
| 20 |
|
|---|
| 21 | run = RawData(data_filename, calib_filename, return_dict = True, do_calibration=True)
|
|---|
| 22 | despike = DRSSpikes()
|
|---|
| 23 | smooth = SlidingAverage(7)
|
|---|
| 24 | cfd = CFD()
|
|---|
| 25 |
|
|---|
| 26 | thr = 3
|
|---|
| 27 | filter_delay = 3
|
|---|
| 28 | search_window_size = 12
|
|---|
| 29 | # shortcut
|
|---|
| 30 | sws = search_window_size
|
|---|
| 31 |
|
|---|
| 32 | plt.ion()
|
|---|
| 33 | fig = plt.figure()
|
|---|
| 34 | fig.hold(True)
|
|---|
| 35 |
|
|---|
| 36 |
|
|---|
| 37 |
|
|---|
| 38 | for event in run:
|
|---|
| 39 | print event['event_id'].value
|
|---|
| 40 | data = event['acal_data']
|
|---|
| 41 | data = despike(data)
|
|---|
| 42 | data_orig = data.copy()
|
|---|
| 43 | data = smooth(data)
|
|---|
| 44 | filtered = cfd(data)
|
|---|
| 45 | filtered = smooth(filtered)
|
|---|
| 46 |
|
|---|
| 47 | # this is a loop over all pixel of this event
|
|---|
| 48 |
|
|---|
| 49 | for dat, fil, orig in zip(data, filtered, data_orig):
|
|---|
| 50 | plt.cla()
|
|---|
| 51 | prod = fil[:-1] * fil[1:]
|
|---|
| 52 | cand = np.where( prod <= 0)[0]
|
|---|
| 53 | if len(cand) == 0:
|
|---|
| 54 | continue
|
|---|
| 55 | # zero crossing with rising edge
|
|---|
| 56 | cross = cand[np.where(fil[cand] < 0)[0]]
|
|---|
| 57 | if len(cross) == 0:
|
|---|
| 58 | continue
|
|---|
| 59 | over_thr = cross[np.where(dat[cross-4] > thr)[0]]
|
|---|
| 60 |
|
|---|
| 61 | # Now since we have these values, we will throw away all those,
|
|---|
| 62 | # which are probably on a falling edge of its predecessor
|
|---|
| 63 | dover = np.diff(over_thr)
|
|---|
| 64 | if len(dover) == 0:
|
|---|
| 65 | good = over_thr
|
|---|
| 66 | else:
|
|---|
| 67 | good = []
|
|---|
| 68 | good.append(over_thr[0])
|
|---|
| 69 | for i in range(len(dover)):
|
|---|
| 70 | if dover[-i-1] > 100:
|
|---|
| 71 | good.append(over_thr[-i-1])
|
|---|
| 72 |
|
|---|
| 73 | good = np.array(sorted(good))
|
|---|
| 74 | # these positions, we just found, do not exactly point to the maximum
|
|---|
| 75 | # of a peak, but the peak will be on the left side of it.
|
|---|
| 76 | # we use the smoothed data to find the position of the local maximum
|
|---|
| 77 | # and then stopre this position and the value of both
|
|---|
| 78 | # the smoothed data and the original data.
|
|---|
| 79 |
|
|---|
| 80 | max_pos = np.zeros( good.shape )
|
|---|
| 81 | max_smoothed = np.zeros( good.shape )
|
|---|
| 82 | max_orig = np.zeros( good.shape )
|
|---|
| 83 |
|
|---|
| 84 | for i in range(len(good)):
|
|---|
| 85 | # We search for a local maximum in a window of size 12
|
|---|
| 86 | if len(dat[good[i]-sws:good[i]]) > 0:
|
|---|
| 87 |
|
|---|
| 88 | max_pos[i] = good[i]-sws + np.argmax(dat[good[i]-sws:good[i]])
|
|---|
| 89 | max_smoothed[i] = dat[max_pos[i]]
|
|---|
| 90 | max_orig[i] = orig[max_pos[i]-filter_delay]
|
|---|
| 91 |
|
|---|
| 92 | plt.plot(max_pos, max_smoothed, 'ro')
|
|---|
| 93 | plt.plot(max_pos, max_orig, 'bo')
|
|---|
| 94 | plt.plot(np.arange(len(dat)), dat, 'k:')
|
|---|
| 95 |
|
|---|
| 96 | raw_input('bla')
|
|---|
| 97 |
|
|---|
| 98 |
|
|---|
| 99 | #output = open(out_filename, 'wb')
|
|---|
| 100 | #cPickle.dump(data_filename, output, -1)
|
|---|
| 101 | #cPickle.dump(calib_filename, output, -1)
|
|---|
| 102 | #cPickle.dump(peak_list, output, -1)
|
|---|
