Index: fact/tools/pyscripts/sandbox/dneise/gain/gain_hdf5.py
===================================================================
--- fact/tools/pyscripts/sandbox/dneise/gain/gain_hdf5.py	(revision 13449)
+++ fact/tools/pyscripts/sandbox/dneise/gain/gain_hdf5.py	(revision 13449)
@@ -0,0 +1,156 @@
+#!/usr/bin/python -tt
+#
+# 
+
+from pyfact     import RawData
+from drs_spikes import DRSSpikes
+from fir_filter import CFD
+from fir_filter import SlidingAverage
+
+import os.path
+import sys
+import numpy as np
+try:
+    import h5py
+except:
+    print 'module h5py needed for this process ... aborting'
+    sys.exit(1)
+
+data_filename = 'data/20111017_010.fits.gz'
+calib_filename = 'data/20111017_006.drs.fits.gz'
+
+run = RawData(data_filename, calib_filename, return_dict = True, do_calibration=True)
+despike = DRSSpikes()
+smooth = SlidingAverage(7)
+cfd = CFD()
+
+thr = 3
+filter_delay = 3
+search_window_size = 12
+# shortcut
+sws = search_window_size
+
+
+# we try to determine the:
+#   * Position
+#   * and Height (in two versions: filtered, and unfiltered data)
+# of dark count peaks in:
+#  a) all events
+#  b) all pixel
+#
+# we limit ourselfs to 10 peaks per pixel max!
+# and 1000 events max!
+#
+# So in order to store this stuff, we prepare an n-dim-array
+
+PEAKS_MAX = 30
+
+# In order not to allocate so much mem, we use hdf5 files for storing this information
+f = h5py.File('20111017_010-006.h5', 'w')
+
+# we write some meta data into the file, so we know afterwards, where this stuff 
+# came from
+f.attrs['data_file'] = os.path.realpath(data_filename)
+f.attrs['calib_file'] = os.path.realpath(calib_filename)
+# this script
+f.attrs['source_python_script'] = open(sys.argv[0], 'r').read()
+# SVN revision, should go here ... but of what? this file??
+f.attrs['settings'] = [
+    'threshold = 3',
+    'SlidingAverage = 7 --> Filter delay = 3'
+    'Search Window Size = 12'
+    'Maximum Number of Peaks, in data sets = 30'
+    'Default value of peak_positions = -1'
+    'Default value of peak_height_unfiltered = -np.inf'
+    'Default value of peak_height_smoothed = -np.inf'
+]
+
+result_peak_positions         = f.create_dataset('peak_positions', 
+                                (run.nevents, run.npix,  PEAKS_MAX), 
+                                dtype=np.int16, 
+                                fillvalue=-1, 
+                                compression='gzip')
+result_peak_unfiltered_height = f.create_dataset('peak_height_unfiltered', 
+                                (run.nevents, run.npix,  PEAKS_MAX), 
+                                dtype=np.float32, 
+                                fillvalue=-np.inf, 
+                                compression='gzip')
+result_peak_smoothed_height   = f.create_dataset('peak_height_smoothed', 
+                                (run.nevents, run.npix,  PEAKS_MAX), 
+                                dtype=np.float32, 
+                                fillvalue=-np.inf, 
+                                compression='gzip')
+
+try:
+
+    for event in run:
+        event_id = event['event_id'].value - 1
+        
+        data = event['acal_data']
+        data = despike(data)
+        data_orig = data.copy()
+        data = smooth(data)
+        filtered = cfd(data)
+        filtered = smooth(filtered)
+
+        # this is a loop over all pixel of this event
+        pixel_id = -1
+        for dat, fil, orig in zip(data, filtered, data_orig):
+            pixel_id += 1
+            print event_id, pixel_id
+
+            prod = fil[:-1] * fil[1:]
+            cand = np.where( prod <= 0)[0]
+            if len(cand) == 0:
+                continue
+            # zero crossing with rising edge
+            cross = cand[np.where(fil[cand] < 0)[0]]
+            if len(cross) == 0:
+                continue
+            over_thr = cross[np.where(dat[cross-4] > thr)[0]]
+
+            # Now since we have these values, we will throw away all those,
+            # which are probably on a falling edge of its predecessor
+            dover = np.diff(over_thr)
+            if len(dover) == 0:
+                good = over_thr
+            else:
+                good = []
+                good.append(over_thr[0])
+                for i in range(len(dover)):
+                    if dover[-i-1] > 100:
+                        good.append(over_thr[-i-1])
+                
+                good = np.array(sorted(good))
+            # 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( good.shape )
+            max_smoothed = np.zeros( good.shape )
+            max_orig     = np.zeros( good.shape )
+        
+            for i in range(len(good)):
+                # We search for a local maximum in a window of size 12
+                if len(dat[good[i]-sws:good[i]]) > 0:
+                
+                    max_pos[i]       = good[i]-sws + np.argmax(dat[good[i]-sws:good[i]])
+                    max_smoothed[i]  = dat[max_pos[i]]
+                    max_orig[i]      = orig[max_pos[i]-filter_delay]
+            
+            result_peak_positions[event_id,pixel_id, :len(max_pos)]         = max_pos
+            result_peak_unfiltered_height[event_id,pixel_id, :len(max_pos)] =max_orig
+            result_peak_smoothed_height[event_id,pixel_id, :len(max_pos)]   = max_smoothed
+
+
+f.close()
+
+except:
+    print 'An Exception occured .. closing file'
+    f.close()
+    raise
+
+
+