Changeset 13439 for fact/tools

Timestamp:

04/24/12 16:52:43 (13 years ago)

Author:

neise

Message:

further steps

File:

• fact/tools/pyscripts/sandbox/dneise/gain/gain.py (modified) (3 diffs)

fact/tools/pyscripts/sandbox/dneise/gain/gain.py

@@ -6 +6 @@
 from drs_spikes import DRSSpikes        
 from fir_filter import CFD
-from extractor import ZeroXing
 from fir_filter import SlidingAverage
 
@@ -16 +15 @@
 import numpy as np
 
-data_filename = 'data/20111017_010.fits.gz'
-calib_filename = 'data/20111017_007.drs.fits.gz'
+data_filename = 'data/20111017_009.fits.gz'
+calib_filename = 'data/20111017_006.drs.fits.gz'
 out_filename = 'test.pkl'
 
 run = RawData(data_filename, calib_filename, return_dict = True, do_calibration=True)
 despike = DRSSpikes()
-smooth = SlidingAverage(8)
+smooth = SlidingAverage(7)
 cfd = CFD()
-find = ZeroXing()
 
-#plt.ion()
+thr = 3
+filter_delay = 3
 
-peak_list = []
-        
-p = Plotter('data[0]')
+plt.ion()
+fig = plt.figure()
+fig.hold(True)
 
-
-#fig = plt.figure()
-#plt.grid(True)
-
-stupid_list = []
 
 for event in run:
@@ -42 +36 @@
     data = event['acal_data']
     data = despike(data)
+    data_orig = data.copy()
     data = smooth(data)
     filtered = cfd(data)
     filtered = smooth(filtered)
-    peak_positions = find(data)
+
+    for dat, fil, orig in zip(data, filtered, data_orig):
+        plt.cla()
+        prod = fil[:-1] * fil[1:]
+        cand = np.where( prod <= 0)[0]
+        # zero crossing with rising edge
+        cross = cand[np.where(fil[cand] < 0)[0]]
+        
+        over_thr = cross[np.where(dat[cross-4] > thr)[0]]
+
+        # Now we have these values, we will throw away all those,
+        # which are probably on a falling edge of its predecessor
+        
+        
+        over = []
+        dover = np.diff(over_thr)
+        
+        for i in range(len(over_thr)):
+            if dover[i] > 100:
+                
+                
+        over_thr = np.array(over)
+        # these positions, we just found, do not exactly point to the maximum
+        # of a peak, but the peak will be on the left side of it.
+        # we use the smoothed data to find the position of the local maximum
+        # and then stopre this position and the value of both
+        # the smoothed data and the original data.
+        
+        max_pos      = np.zeros( over_thr.shape )
+        max_smoothed = np.zeros( over_thr.shape )
+        max_orig     = np.zeros( over_thr.shape )
     
-    # find peak heights
-    # and get rid of too small peaks
-#    good_peaks = []
-#    for idx,pixel in enumerate(peak_positions):
-#        good_peaks.append([])
-#        for peak_pos in pixel:
-#            peak_height = data[idx][int(peak_pos)]
-#            stupid_list.append(peak_height)
-#            if peak_height > 4.0:
-#                good_peaks[-1].append((peak_pos,peak_height)) 
-#
-#    plt.hist(np.array(stupid_list))
+        for i in range(len(over_thr)):
+            # We search for a local maximum in a window of size 12
+            if len(dat[over_thr[i]-12:over_thr[i]]) > 0:
+            
+                max_pos[i]       = over_thr[i]-12 + np.argmax(dat[over_thr[i]-12:over_thr[i]])
+                max_smoothed[i]  = dat[max_pos[i]]
+                max_orig[i]       = orig[max_pos[i]-filter_delay]
 
-    p((data[0],filtered[0]))
-    print 'peak_positions:', map(int, peak_positions[0])
-    
-    peak_pussies = map(int, peak_positions[0])
-    peak_h = []
-    for i in peak_pussies:
-        peak_h.append(data[0][i-8])
-
-
-    pp = Plotter('peaks', x=peak_pussies)
-    pp(peak_h)
-    
-
-    ret = raw_input('quit?')
-    if ret=='q':
-        break
-
-#    peak_list.append( good_peaks )
+        plt.plot(max_pos, max_smoothed, 'ro')
+        plt.plot(max_pos, max_orig, 'bo')
+        plt.plot(np.arange(len(dat)), dat, 'k:')
+        
+        raw_input('bla')