| 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 extractor import ZeroXing
|
|---|
| 9 | from fir_filter import SlidingAverage
|
|---|
| 10 |
|
|---|
| 11 | from plotters import Plotter
|
|---|
| 12 |
|
|---|
| 13 | import sys
|
|---|
| 14 | import cPickle
|
|---|
| 15 | import matplotlib.pyplot as plt
|
|---|
| 16 | import numpy as np
|
|---|
| 17 |
|
|---|
| 18 | data_filename = 'data/20111017_009.fits.gz
|
|---|
| 19 | calib_filename = 'data/20111017_006.drs.fits.gz'
|
|---|
| 20 | out_filename = 'test.pkl'
|
|---|
| 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 | find = ZeroXing()
|
|---|
| 27 |
|
|---|
| 28 | plt.ion()
|
|---|
| 29 |
|
|---|
| 30 | peak_list = []
|
|---|
| 31 |
|
|---|
| 32 | #p = Plotter('data[0]')
|
|---|
| 33 |
|
|---|
| 34 |
|
|---|
| 35 | #fig = plt.figure()
|
|---|
| 36 | #plt.grid(True)
|
|---|
| 37 |
|
|---|
| 38 | stupid_list = []
|
|---|
| 39 | thr = 3
|
|---|
| 40 |
|
|---|
| 41 | fig = plt.figure()
|
|---|
| 42 | fig.hold(True)
|
|---|
| 43 |
|
|---|
| 44 | for event in run:
|
|---|
| 45 | print event['event_id'].value
|
|---|
| 46 | data = event['acal_data']
|
|---|
| 47 | origorig = data.copy()
|
|---|
| 48 | data = despike(data)
|
|---|
| 49 | data_orig = data.copy()
|
|---|
| 50 | data = smooth(data)
|
|---|
| 51 | filtered = cfd(data)
|
|---|
| 52 | filtered = smooth(filtered)
|
|---|
| 53 |
|
|---|
| 54 | for dat, fil, orig, zzz in zip(data, filtered, data_orig, origorig):
|
|---|
| 55 | plt.cla()
|
|---|
| 56 | prod = fil[:-1] * fil[1:]
|
|---|
| 57 | cand = np.where( prod <= 0)[0]
|
|---|
| 58 | # zero crossing with rising edge
|
|---|
| 59 | cross = cand[np.where(fil[cand] < 0)[0]]
|
|---|
| 60 |
|
|---|
| 61 | over_thr = cross[np.where(dat[cross-4] > thr)[0]]
|
|---|
| 62 |
|
|---|
| 63 | maxima = np.zeros( over_thr.shape )
|
|---|
| 64 | maxima_pos =np.zeros( over_thr.shape )
|
|---|
| 65 | omaxima = np.zeros( over_thr.shape )
|
|---|
| 66 | omaxima_pos =np.zeros( over_thr.shape )
|
|---|
| 67 |
|
|---|
| 68 | for i in range(len(over_thr)):
|
|---|
| 69 | if len(dat[over_thr[i]-12:over_thr[i]]) > 0:
|
|---|
| 70 | maxima_pos[i] = over_thr[i]-12 + np.argmax(dat[over_thr[i]-12:over_thr[i]])
|
|---|
| 71 | maxima[i] = dat[maxima_pos[i]]
|
|---|
| 72 | omaxima_pos[i] = over_thr[i]-12 + np.argmax(orig[over_thr[i]-12:over_thr[i]])
|
|---|
| 73 | omaxima[i] = orig[maxima_pos[i]-3]
|
|---|
| 74 |
|
|---|
| 75 | x_values = np.array(range(len(dat)))
|
|---|
| 76 |
|
|---|
| 77 | plt.plot( x_values-3, dat , 'k:')
|
|---|
| 78 | plt.plot( x_values, orig , 'y:')
|
|---|
| 79 | plt.plot( x_values, zzz , 'r.')
|
|---|
| 80 |
|
|---|
| 81 | plt.plot( maxima_pos-3, maxima, 'bo')
|
|---|
| 82 | plt.plot( maxima_pos-3, omaxima, 'yo')
|
|---|
| 83 | plt.plot( over_thr-3, dat[over_thr], 'ro')
|
|---|
| 84 |
|
|---|
| 85 |
|
|---|
| 86 |
|
|---|
| 87 | ret = raw_input('quit?')
|
|---|
| 88 | if ret=='q':
|
|---|
| 89 | break
|
|---|
| 90 |
|
|---|
| 91 | # peak_list.append( good_peaks )
|
|---|
| 92 |
|
|---|
| 93 |
|
|---|
| 94 | #output = open(out_filename, 'wb')
|
|---|
| 95 | #cPickle.dump(data_filename, output, -1)
|
|---|
| 96 | #cPickle.dump(calib_filename, output, -1)
|
|---|
| 97 | #cPickle.dump(peak_list, output, -1)
|
|---|
