Index: fact/tools/pyscripts/sandbox/dneise/spike_studies/drs_spikes.py
===================================================================
--- fact/tools/pyscripts/sandbox/dneise/spike_studies/drs_spikes.py	(revision 13632)
+++ fact/tools/pyscripts/sandbox/dneise/spike_studies/drs_spikes.py	(revision 13632)
@@ -0,0 +1,120 @@
+#!/usr/bin/python -tt
+#
+# Werner Lustermann
+# ETH Zurich
+#
+import numpy as np
+
+import fir_filter as fir
+
+class DRSSpikes(object):
+    """ remove spikes (single or double false readings) from DRS4 data
+    Strategy:
+    * filter the data, removing the signal, thus spike(s) are clearly visible
+    * search single and double spikes
+    * replace the spike by a value derived from the neighbors
+    
+    """
+    
+    def __init__(self, threshold=7., 
+                 single_pattern=np.array( [-0.5, 1.0, -0.5]) ,
+                 double_pattern=np.array([-1., 1., 1., -1.]), 
+                 user_action=lambda DRSSpikes_instance: None,
+                 debug = False):
+        """ initialize spike filter 
+        template_single: template of a single slice spike
+        template_double: template of a two slice spike
+
+        """
+
+        self.threshold = threshold
+        self.single_pattern = single_pattern * threshold
+        self.double_pattern = double_pattern * threshold
+        
+        self.remove_signal = fir.RemoveSignal()
+
+        self.user_action = user_action
+        self.debug = debug
+
+    def __call__(self, data):
+
+        self.row, self.col = data.shape
+        indicator = self.remove_signal(data)
+        a = indicator.flatten()
+        singles = []
+        doubles = []
+
+        # a spike in the first or last channel is considered as a filter artefact
+        candidates = np.where(a[1:-2] > self.threshold)
+        # candidates = np.where(a[1:1022] > self.threshold)
+        cc = candidates[0]
+        #print 'cc: ', cc
+        #: find single spikes
+        p = self.single_pattern * np.sign( self.single_pattern )
+        for i, can in enumerate( zip(a[cc], a[cc+1], a[cc+2]) ):
+            #print 'can : p', can, p
+            can = can * np.sign(self.single_pattern)
+            if np.all(can > p):
+                singles.append(cc[i])
+
+        #: find double spikes
+        p = self.double_pattern * np.sign( self.double_pattern )
+        for i, can in enumerate( zip(a[cc], a[cc+1], a[cc+2], a[cc+3]) ):
+            #print 'data:  ', [data[0,cc[i]+k] for k in range(3)]
+            #print 'can : p', can, p
+            can = can * np.sign(self.double_pattern)
+            if np.all(can > p):
+                doubles.append(cc[i])
+
+        if self.debug:
+            print 'singles: ', singles
+            print 'doubles: ', doubles
+        
+        # make relevant information a member,
+        self.singles = singles
+        self.doubles = doubles 
+        self.indicator = a
+        self.candidates = cc
+
+        self.user_action(self)
+
+        data = self.remove_single_spikes(singles, data)
+        data = self.remove_double_spikes(doubles, data)
+        return data
+
+    def remove_single_spikes(self, singles, data):
+        data = data.flatten() 
+        for spike in singles:
+            data[spike] = (data[spike-1] + data[spike+1]) / 2.
+        return data.reshape(self.row, self.col)
+    
+    def remove_double_spikes(self, doubles, data):
+        data = data.flatten() 
+        for spike in doubles:
+            data[spike:spike+2] = (data[spike-1] + data[spike+2]) / 2.
+        return data.reshape(self.row, self.col)
+
+
+def _test():
+  
+    a = np.ones((3,12)) * 3.
+    a[0,3] = 7.
+    a[1,7] = 14.
+    a[1,8] = 14.
+    a[2,4] = 50.
+    a[2,5] = 45.
+    a[2,8] = 20.
+    
+    print a 
+
+    SpikeRemover = DRSSpikes(3., debug=True)
+    print 'single spike pattern ', SpikeRemover.single_pattern
+    print 'double spike pattern ', SpikeRemover.double_pattern
+    afilt = SpikeRemover(a)
+    print afilt
+    
+if __name__ == '__main__':
+    """ test the class """
+    _test()
+    
+    
Index: fact/tools/pyscripts/sandbox/dneise/spike_studies/drs_spikes_2d.py
===================================================================
--- fact/tools/pyscripts/sandbox/dneise/spike_studies/drs_spikes_2d.py	(revision 13632)
+++ fact/tools/pyscripts/sandbox/dneise/spike_studies/drs_spikes_2d.py	(revision 13632)
@@ -0,0 +1,120 @@
+class DRSSpikes_2D(object):
+    """ remove spikes (single or double false readings) from DRS4 data
+    Strategy:
+    * filter the data, removing the signal, thus spike(s) are clearly visible
+    * search single and double spikes
+    * replace the spike by a value derived from the neighbors
+    
+    """
+    
+    def __init__(self, threshold=7., 
+                 single_pattern=np.array( [-0.5, 1.0, -0.5]) ,
+                 double_pattern=np.array([-1., 1., 1., -1.]), 
+                 user_action=lambda self, data: None,
+                 debug = False):
+        """ initialize spike filter 
+        template_single: template of a single slice spike
+        template_double: template of a two slice spike
+        """
+
+        self.threshold = threshold
+        self.single_pattern = single_pattern * threshold
+        self.double_pattern = double_pattern * threshold
+        
+        self.remove_signal = fir.RemoveSignal()
+
+        self.user_action = user_action
+        self.debug = debug
+
+    def __call__(self, data):
+        # shortcuts
+        row, col = data.shape
+        thr = self.threshold
+        
+        # these are: lists, which will contain positiones of spikes
+        # lets see if this is feasible
+        self.singles = []
+        self.doubles = []
+        singles = self.singles
+        doubles = self.doubles
+
+        # indi means indicator, i.e. a filter output, which indicates, where spikes
+        # are positioned in the unfiltered data
+        # indi is delayed w.r.t. data by 1
+        self.indi = self.remove_signal(data)
+        indi = self.indi
+
+        # cand (candidates), is a tuple of two equal length np.arrays
+        # each pair of array elements can be understood as coordinates, pointing out
+        # where the condition was fullfilled in indi. 
+        # e.g. in pixel = cand[0][0] around slice = cand[1][0]
+        # there is probably a spike
+        #
+        # for double spikes, two neighboring slices fulfill the condition
+        # which lead to something like: 
+        # cand =( array([ ... 3, 3, ... ]) , array([ ... 102, 103 ...]) )
+        cand = np.where(indi[:,1:-2] > thr)
+        self.cand = cand
+        # in order to verify, that the candidate is really a single or double
+        # spike, we compare the spike with a 3 or 4 slices pattern.
+        # therefor we want to slice out 4 slices out of indi, where ever the 
+        # condition was fullfilled
+        #
+        # note: since indi was sliced in the np.where statement, 
+        # the resulting cand coordinates are reduced by 1 in the slice coordinate
+        # this is actually what we want, since a spike has a distinctive low-high-low 
+        # pattern in the indicator. So we *want* the indicator slices to be shifted 1
+        # to the left
+        # and in addition, by pure chance, the coordinates in cand[1] point directly 
+        # to the real spike in data, since indi was delayed by one anyway.
+        cand_slices = np.empty( (len(cand[0]), 4), dtype=np.int )
+        self.cand_slices = cand_slices
+        for i in range(4):
+            cand_slices[i] = indi[ (cand[0], cand[1]+i )
+            
+        # search for single spikes
+        sp = self.single_pattern * np.sign( self.single_pattern )
+        for i, can in enumerate(cand_slices[:-1]):
+            can *= np.sign(sp)
+            if np.all( can > sp):
+                singles.append( (cand[0][i],can[1][i]) )
+
+        # I guess in principle it is possible, that a candidate looks like a 
+        # single and like a double, ... nut with the current patters
+        # but in case one changes the patterns ... then it might happen.
+        # In addition the treatment of double spikes is maybed not smart:
+        #   In case both parts of a double spike fulfill the 1st condition
+        #       only the first candidate will fulfill the 2nd condition
+        #   In case only the first part fulfilled the 1st conition, then
+        #       we are fine
+        #   In case only the second part triggered the first time, then
+        #       we sliced out the wrong piece and it wouldn't fulfull the 2nd anyway.
+        # 
+        #   This means, in case there are neighboring hits in cand,
+        #   The 2nd neighbor will never fulfill the 2nd condition.
+
+        # search for double spikes
+        dp = self.double_pattern * np.sign( self.double_pattern )
+        for i, can in enumerate( cand_slices ):
+            can *= np.sign(dp)
+            if np.all(can > dp):
+                doubles.append( (cand[0][i],can[1][i]) )
+
+        self.user_action(self, data)
+
+        data = self.remove_single_spikes(singles, data)
+        data = self.remove_double_spikes(doubles, data)
+        return data
+
+    def remove_single_spikes(self, singles, data):
+        for spike in singles:
+            data[spike[0],spike[1]] = (data[spike[0],spike[1]-1] + data[spike[0],spike[1]+1]) / 2.
+        return data
+    
+    def remove_double_spikes(self, doubles, data):
+        for spike in doubles:
+            data[spike[0],spike[1]:spike[1]+2] = (data[spike[0],spike[1]-1] + data[spike[0],spike[1]+2]) / 2.
+        return data
+
+
+
Index: fact/tools/pyscripts/sandbox/dneise/spike_studies/save_spikes.py
===================================================================
--- fact/tools/pyscripts/sandbox/dneise/spike_studies/save_spikes.py	(revision 13632)
+++ fact/tools/pyscripts/sandbox/dneise/spike_studies/save_spikes.py	(revision 13632)
@@ -0,0 +1,41 @@
+#!/usr/bin/python -tti
+
+# script, which saves a lot of information about 
+# DRS spikes in a file, for further analysis.
+from pyfact import RawData
+from drs_spikes import DRSSpikes
+import numpy as np
+
+# I need a callback function for the call of DRSSpikes
+# this function should be a member of RawData, 
+# so that it knows all necessary information
+# In order to do this, I should inherit from RawData
+# and just add one single function.
+class RD_spike_ana( RawData ):
+    def callback(self, drsspikes ):
+        print 'Event ID', self.event_id
+        print 'singles:'
+        for s in drsspikes.singles:
+            print np.unravel_index( s, (drsspikes.row, drsspikes.col) )
+
+def spike_callback(drsspikes):
+    print drsspikes.__dict__
+
+
+data_file_name = '/home/dominik/20120223_205.fits.gz'
+calib_file_name = '/home/dominik/20120223_202.drs.fits.gz'
+outfile = None
+
+run = RD_spike_ana( data_file_name, calib_file_name, return_dict=True)
+despike = DRSSpikes(user_action = run.callback)
+
+
+
+
+for event in run:
+    data  = event['acal_data'][:,12:-12] 
+    ddata = despike(data)
+
+    ret = raw_input('q for quit')
+    if 'q' in ret:
+        break
Index: fact/tools/pyscripts/sandbox/dneise/spike_studies/spike_ana.py
===================================================================
--- fact/tools/pyscripts/sandbox/dneise/spike_studies/spike_ana.py	(revision 13632)
+++ fact/tools/pyscripts/sandbox/dneise/spike_studies/spike_ana.py	(revision 13632)
@@ -0,0 +1,43 @@
+#!/usr/bin/python -tti
+
+# script, which saves a lot of information about 
+# DRS spikes in a file, for further analysis.
+from pyfact import RawData
+from drs_spikes import DRSSpikes
+import numpy as np
+from plotters import Plotter
+import scipy.signal as signal
+# I need a callback function for the call of DRSSpikes
+# this function should be a member of RawData, 
+# so that it knows all necessary information
+# In order to do this, I should inherit from RawData
+# and just add one single function.
+data_file_name = '/home/dominik/20120223_205.fits.gz'
+calib_file_name = '/home/dominik/20120223_202.drs.fits.gz'
+outfile = None
+
+run = RawData( data_file_name, calib_file_name, return_dict=True)
+despike = DRSSpikes()
+
+single_pattern = np.array( (-1,2,-1) )
+double_pattern = np.array( (-1,1,1,-1) )
+
+bins = None
+
+histplot = Plotter('hist')
+cplot = Plotter('corr')
+for event in run:
+    print event['event_id']
+    data  = event['acal_data'][:,12:-12]
+    
+    for pixel,sig in enumerate(data):
+        corr = signal.correlate( sig, single_pattern, 'same')
+        #cplot.name = corr+'_'+str(event['event_id']+'_'+str(pixel)
+        cplot(corr)
+        hist = np.histogram( corr, bins=1000, range=(0,2000) )
+        if bins == None:
+            bins = hist[0]
+        else:
+            bins += hist[0]
+            
+    histplot(bins)
Index: fact/tools/pyscripts/sandbox/dneise/spike_studies/zerosearch.py
===================================================================
--- fact/tools/pyscripts/sandbox/dneise/spike_studies/zerosearch.py	(revision 13632)
+++ fact/tools/pyscripts/sandbox/dneise/spike_studies/zerosearch.py	(revision 13632)
@@ -0,0 +1,224 @@
+#!/usr/bin/python -tt
+#
+# Dominik Neise, Werner Lustermann
+# TU Dortmund, ETH Zurich
+#
+import numpy as np
+
+class ZeroSearch( object ):
+    """ class closely related to the methods in rootmacros/zerosearch.C and .h
+    """
+    def __init__(self, slope=1, name = 'ZeroSearch'):
+        
+    def __call__(self, data):
+        """ returns list of lists """
+        
+        for pid, pixel in enumerate(data):
+            
+            # candidates for zero crossings are, where two neighboring 
+            # cells have opposite signs, i.e. their product is < 0
+            product = pixel[:-1] * pixel[1:]
+            candidates = np.where(product < 0)[0]
+            
+            if slope>0:
+                candidates = candidates[np.where(pixel[candidates] > 0)[0]]
+            else:
+                candidates = candidates[np.where(pixel[candidates] < 0)[0]]
+                
+
+class GlobalMaxFinder(object):
+    """ Pulse Extractor
+        Finds the global maximum in the given window.
+        (Best used with filtered data)
+    """
+    
+    def __init__(self, min=30, max=250 , name = 'GlobalMaxFinder'):
+        """ initialize search Window
+        
+        """
+        self.__module__="extractor"
+        self.min = min
+        self.max = max
+        self.name = name
+        
+    def __call__(self, data):
+        if data.ndim > 1:
+            time = np.argmax( data[ : , self.min:self.max ], 1)
+            amplitude = np.max( data[ : , self.min:self.max], 1)
+        else:
+            time = np.argmax( data[self.min:self.max])
+            amplitude = np.max( data[self.min:self.max])
+        return amplitude, time+self.min
+
+    def __str__(self):
+        s = self.name + '\n'
+        s += 'window:\n'
+        s += '(min,max) = (' + str(self.min) + ',' + str(self.max) + ')'
+        return s
+        
+    def test(self):
+        pass
+
+
+class WindowIntegrator(object):
+    """ Integrates in a given intergration window around the given position
+    """
+    
+    def __init__(self, min=13, max=23 , name = 'WindowIntegrator'):
+        """ initialize integration Window
+        """
+        self.__module__="extractor"
+        self.min = min
+        self.max = max
+        self.name = name
+        
+    def __call__(self, data, pos):
+        integral = np.empty( data.shape[0] )
+        for pixel in range( data.shape[0] ):
+            integral[pixel] = data[pixel, (pos[pixel]-self.min):(pos[pixel]+self.max)].sum()
+        return integral
+
+    def __str__(self):
+        s = self.name + '\n'
+        s += 'window:\n'
+        s += '(min,max) = (' + str(self.min) + ',' + str(self.max) + ')'
+        return s    
+
+class FixedWindowIntegrator(object):
+    """ Integrates in a given intergration window
+    """
+    
+    def __init__(self, min=55, max=105 , name = 'FixedWindowIntegrator'):
+        """ initialize integration Window
+        """
+        self.__module__="extractor"
+        self.min = min
+        self.max = max
+        self.name = name
+        
+    def __call__(self, data):
+        integral = np.empty( data.shape[0] )
+        for pixel in range( data.shape[0] ):
+            integral[pixel] = data[pixel, self.min:self.max].sum()
+        return integral
+
+    def __str__(self):
+        s = self.name + '\n'
+        s += 'window:\n'
+        s += '(min,max) = (' + str(self.min) + ',' + str(self.max) + ')'
+        return s    
+
+class ZeroXing(object):
+    """ Finds zero crossings in given data
+        (should be used on CFD output for peak finding)
+        returns list of lists of time_of_zero_crossing
+    """
+    def __init__(self, slope=1, name = 'ZeroXing'):
+        self.__module__="extractor"
+        if (slope >= 0):
+            self.slope = 1  # search for rising edge crossing
+        elif (slope < 0):
+            self.slope = -1 # search for falling edge crossing
+        self.name = name
+
+
+    def __call__(self, data, zero_level = 0):
+        all_hits = []
+        for pix_data in data:
+            hits = []
+            for i in range( data.shape[1]-1 ):
+                dat = pix_data[i] - zero_level
+                next_dat = pix_data[i+1] - zero_level
+                if ( self.slope > 0 ):
+                    if ( dat > 0 ):
+                        continue
+                else:
+                    if ( dat < 0):
+                        continue
+                if ( dat * next_dat <= 0 ):
+                    # interpolate time of zero crossing with 
+                    # linear polynomial: y = ax + b
+                    a = (next_dat - dat) / ((i+1) - i)
+                    time = -1.0/a * dat + i
+                    hits.append(time)
+            all_hits.append(hits)
+        return all_hits
+
+    def __str__(self):
+        s = self.name + '\n'
+        if (self.slope == 1):
+            s += 'search for rising edge crossing.\n'
+        else:
+            s += 'search for falling edge crossing.\n'
+        return s
+
+
+
+def _test_GlobalMaxFinder():
+    gmf = GlobalMaxFinder(30,250)
+    print gmf
+    amplitude, time = gmf(event)
+    if abs(amplitude.mean() - 10) < 0.5:
+        print "Test 1: OK GlobalMaxFinder found amplitude correctly", amplitude.mean()
+    if abs(time.mean() - 65) < 2:
+        print "Test 1: OK GlobalMaxFinder found time correctly", time.mean()
+    else:
+        print "BAD: time mean:", time.mean()
+
+def _test_FixedWindowIntegrator():
+    fwi = FixedWindowIntegrator(50,200)
+    print fwi
+    integral = fwi(event)
+    #value of integral should be: 150*bsl + 8*10/2 + 100*10/2 = 465
+    if abs( integral.mean() - 465) < 2:
+        print "Test 2: OK FixedWindowIntegrator found integral correctly", integral.mean()
+    else:
+        print "Test 2:  X FixedWindowIntegrator integral.mean failed:", integral.mean()
+
+def _test_ZeroXing():
+    cfd = CFD()
+    sa = SlidingAverage(8)
+    print sa
+    cfd_out = sa(event)
+    cfd_out = cfd(cfd_out   )
+    cfd_out = sa(cfd_out)
+    zx = ZeroXing()
+    print zx
+    list_of_list_of_times = zx(cfd_out)
+    times = []
+    for list_of_times in list_of_list_of_times:
+        times.extend(list_of_times)
+    times = np.array(times)
+    
+    hist,bins = np.histogram(times,3000,(0,300))
+    most_probable_time = np.argmax(hist)
+    print 'most probable time of zero-crossing', most_probable_time/10.
+    print 'this includes filter delays ... for average filter setting 8 this turns out to be 78.8 most of the time'
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    """ test the extractors """
+    
+    # Generate a fake event, with a triangular pulse at slice 65
+    sg = SignalGenerator()
+    pulse_str = 'len 300 bsl -0.5 noise 0.5 triangle 65 10 8 100' 
+    pulse = sg(pulse_str)
+    event = []
+    for i in range(1440):
+        event.append(sg(pulse_str))
+    event = np.array(event)
+    print 'test event with 1000 pixel generated, like this:'
+    print pulse_str
+    print
+    
+    print '_test_GlobalMaxFinder()'
+    _test_GlobalMaxFinder()
+    print 
+    print
+    print '_test_FixedWindowIntegrator()'
+    _test_FixedWindowIntegrator()
+    print
+    print
+    print '_test_ZeroXing()'
+    _test_ZeroXing()
+    print