Index: fact/tools/pyscripts/sandbox/dneise/spikes/ana.py
===================================================================
--- fact/tools/pyscripts/sandbox/dneise/spikes/ana.py	(revision 14158)
+++ fact/tools/pyscripts/sandbox/dneise/spikes/ana.py	(revision 14171)
@@ -5,4 +5,5 @@
 import sys
 import numpy as np
+from ctypes import *
 import os
 import os.path as path
@@ -12,94 +13,29 @@
 dfn = sys.argv[1]
 cfn = sys.argv[2]
+run = RawData( dfn, cfn) 
+
 
 file_base_name = path.splitext(path.basename(dfn))[0]
-root_filename = '/home_nfs/isdc/neise/' + file_base_name + time.strftime('%Y%m%d_%H%M%S') +'.root'
-
-run = RawData( dfn, cfn) 
+root_filename = '/home_nfs/isdc/neise/' + file_base_name + '_' + time.strftime('%Y%m%d_%H%M%S') +'_spikeana.root'
 
 rootfile = TFile(root_filename, "RECREATE")
 
-canv = TCanvas('canv', 'test canvas',10,10,400,400)
-canv.Divide(2,4)
+baum = TTree('spiketree', 'spike ana tree')
+
+# prepare some vars for the tree
+chid               =  c_int(0)
+startcell          =  c_int(0)
+number_of_singles  =  c_int(0)
+position_of_spikes_in_logical_pipeline   =  np.zeros(100, np.int32)
+position_of_spikes_in_physical_pipeline  =  np.zeros(100, np.int32)
+time_to_previous   =  c_long(0)
 
 
-
-hs_log_chid = TH2F("hs_log_chid",
-    "singles - logical pipe vs chid",
-    1440, -0.5, 1439.5, 
-    300, -0.5, 299.5)
-hd_log_chid = TH2F("hd_log_chid",
-    "doubles - logical pipe vs chid",
-    1440, -0.5, 1439.5, 
-    300, -0.5, 299.5)
-
-hs_phys_log = TH2F("hs_phys_log",
-    "singles logical pipe vs physical pipe",
-    1024, -0.5, 1023.5, 
-    300, -0.5, 299.5)
-
-hd_phys_log = TH2F("hd_phys_log",
-    "doubles logical pipe vs physical pipe",
-    1024, -0.5, 1023.5, 
-    300, -0.5, 299.5)
-
-hs_phys_sc = TH2F("hs_phys_sc",
-    "singles physical pipe vs startcell",
-    1024, -0.5, 1023.5, 
-    1024, -0.5, 1023.5)
-    
-hd_phys_sc = TH2F("hd_phys_sc",
-    "doubles physical pipe vs startcell",
-    1024, -0.5, 1023.5, 
-    1024, -0.5, 1023.5)
-    
-hs_log_sc = TH2F("hs_log_sc",
-    "singles logical pipe vs startcell",
-    1024, -0.5, 1023.5,
-    300, -0.5, 299.5 
-    )
-    
-hd_log_sc = TH2F("hd_log_sc",
-    "doubles logical pipe vs startcell",
-    1024, -0.5, 1023.5,
-    300, -0.5, 299.5, 
-    )
-
-rootfile.mkdir('per_chid')
-rootfile.cd('per_chid')
-
-hd_log_sc_pp = []
-for chid in range(1440):
-    hd_log_sc_pp.append(
-        TH2F("chid_"+str(chid),
-            "doubles - log pipe vs startcell - chid:" +str(chid),
-            1024, -0.5, 1023.5, 
-            300, -0.5, 299.5)
-        )
-        
-rootfile.cd('..')
-
-    
-singlebaum = TTree('singlebaum', 'spike ana tree')
-doublebaum = TTree('doublebaum', 'spike ana tree')
-
-
-# prepare some vars for the tree
-chid = 0
-startcell = 0
-number_of_singles = 0
-number_of_doubles = 0
-single_positions = np.zeros(10)
-double_positions = np.zeros(10)
-singlebaum.Branch('chid',chid,'chid/I')
-singlebaum.Branch('sc',startcell,'sc/I')
-singlebaum.Branch('n',number_of_singles,'n/I')
-singlebaum.Branch('spikes',single_positions,'spikes/I')
-
-doublebaum.Branch('chid',chid,'chid/I')
-doublebaum.Branch('sc',startcell,'sc/I')
-doublebaum.Branch('n',number_of_doubles,'n/I')
-doublebaum.Branch('spikes',double_positions,'spikes/I')
-
+baum.Branch('chid',chid,'chid/I')
+baum.Branch('sc',startcell,'sc/I')
+baum.Branch('n',number_of_singles,'n/I')
+baum.Branch('logpos',position_of_spikes_in_logical_pipeline,'logpos[n]/I')
+baum.Branch('physpos',position_of_spikes_in_physical_pipeline,'physpos[n]/I')
+baum.Branch('time',time_to_previous,'time/I')
 
 
@@ -116,4 +52,6 @@
             
     
+despike = DRSSpikes(user_action = spikecallback)
+
 def mars_spikes( data ):
     """
@@ -167,5 +105,10 @@
     return singles, doubles        
 
-despike = DRSSpikes(user_action = spikecallback)
+
+
+
+
+event = run.next()
+bt_old = event['board_times'].copy
 
 for event in run:
@@ -173,85 +116,21 @@
     s, d = mars_spikes(data)
     sc = event['start_cells']
+    bt = event['board_times'].copy()
+    
+    if len(s) >0 :
+        chid.value  = s[0][0]
+        startcell.value = sc[ chid.value ]
+        number_of_singles.value = len(s)
+        time.value = bt[s[0][0]/9]-bt_old[s[0][0]/9]
+        for i in s:
+            log = i[1]
+            phys  = (startcell.value+log)%1024
+            position_of_spikes_in_logical_pipeline[i] = log
+            position_of_spikes_in_physical_pipeline[i] = phys
+        baum.Fill()
+    
+    bt_old = event['board_times'].copy()
 
-    if len(s) >0 :
-        for i in s:
-            log   = i[1]
-            chid  = i[0]
-            stace = sc[ chid ]
-            phys  = (stace+log)%1024
-            
-            hs_log_chid.Fill( chid, log)
-            hs_phys_log.Fill( phys, log)
-            hs_phys_sc.Fill( stace, phys)
-            hs_log_sc.Fill( stace, log)
-            
-            
-    if len(d) >0 :
-        for i in d:
-            log   = i[1]
-            chid  = i[0]
-            stace = sc[ chid ]
-            phys  = (stace+log)%1024
+baum.Write()
 
-            hd_log_chid.Fill( chid, log)
-            hd_phys_log.Fill( phys, log)
-            hd_phys_sc.Fill( stace, phys)
-            hd_log_sc.Fill( stace, log)
-            hd_log_sc_pp[chid].Fill(stace, log)
-    
-    
-    
-    if event['event_id'] % 100 == 0:
-        print event['event_id']
-        canv.cd(1)
-        hs_log_chid.Draw("COLZ")
-        canv.cd(2)
-        hd_log_chid.Draw("COLZ")
-        canv.cd(3)
-        hs_phys_log.Draw("COLZ")
-        canv.cd(4)
-        hd_phys_log.Draw("COLZ")
-    
-        canv.cd(5)
-        hs_phys_sc.Draw("COLZ")
-        canv.cd(6)
-        hd_phys_sc.Draw("COLZ")
-        canv.cd(7)
-        hs_log_sc.Draw("COLZ")
-        canv.cd(8)
-        hd_log_sc.Draw("COLZ")
-        canv.Modified()
-        canv.Update()
-        
-        
-
-    if event['event_id'] % 1000 == 0:
-        print event['event_id']
-        hs_log_chid.Write('',TObject.kOverwrite)
-        hd_log_chid.Write('',TObject.kOverwrite)
-        hs_phys_log.Write('',TObject.kOverwrite)
-        hd_phys_log.Write('',TObject.kOverwrite)
-        hs_phys_sc.Write('',TObject.kOverwrite)
-        hd_phys_sc.Write('',TObject.kOverwrite)
-        hs_log_sc.Write('',TObject.kOverwrite)
-        hd_log_sc.Write('',TObject.kOverwrite)
-        
-        rootfile.cd('per_chid')
-        for h in hd_log_sc_pp:
-            h.Write('',TObject.kOverwrite)
-        rootfile.cd('..')
-    
-hs_log_chid.Write('',TObject.kOverwrite)
-hd_log_chid.Write('',TObject.kOverwrite)
-hs_phys_log.Write('',TObject.kOverwrite)
-hd_phys_log.Write('',TObject.kOverwrite)
-hs_phys_sc.Write('',TObject.kOverwrite)
-hd_phys_sc.Write('',TObject.kOverwrite)
-hs_log_sc.Write('',TObject.kOverwrite)
-hd_log_sc.Write('',TObject.kOverwrite)
-        
-rootfile.cd('per_chid')
-for h in hd_log_sc_pp:
-    h.Write()
-rootfile.cd('..')    
 rootfile.Close()
Index: fact/tools/pyscripts/sandbox/dneise/spikes/ana_old.py
===================================================================
--- fact/tools/pyscripts/sandbox/dneise/spikes/ana_old.py	(revision 14171)
+++ fact/tools/pyscripts/sandbox/dneise/spikes/ana_old.py	(revision 14171)
@@ -0,0 +1,279 @@
+#!/usr/bin/python -tti
+
+from pyfact import RawData
+from drs_spikes import DRSSpikes
+import sys
+import numpy as np
+from ctypes import *
+import os
+import os.path as path
+from ROOT import TFile, TCanvas, TH2F, TTree, TStyle, TObject
+import time
+
+dfn = sys.argv[1]
+cfn = sys.argv[2]
+run = RawData( dfn, cfn) 
+
+
+file_base_name = path.splitext(path.basename(dfn))[0]
+root_filename = '/home_nfs/isdc/neise/' + file_base_name + '_' + time.strftime('%Y%m%d_%H%M%S') +'_spikeana.root'
+
+rootfile = TFile(root_filename, "RECREATE")
+
+singlebaum = TTree('singlebaum', 'spike ana tree')
+
+# prepare some vars for the tree
+chid               =  c_int(0)
+startcell          =  c_int(0)
+number_of_singles  =  c_int(0)
+position_of_spikes_in_logical_pipeline   =  np.zeros(100, np.int32)
+position_of_spikes_in_physical_pipeline   =  np.zeros(100, np.int32)
+time_to_previous   =  c_long(0)
+
+singlebaum.Branch('chid',chid,'chid/I')
+singlebaum.Branch('sc',startcell,'sc/I')
+singlebaum.Branch('n',number_of_singles,'n/I')
+singlebaum.Branch('logpos',position_of_spikes_in_logical_pipeline,'logpos[n]/I')
+singlebaum.Branch('physpos',position_of_spikes_in_physical_pipeline,'physpos[n]/I')
+singlebaum.Branch('time',time_to_previous,'time/I')
+
+
+def spikecallback(candidates, singles, doubles, data, ind):
+    if len(singles) >0 :
+        for s in singles:
+            s = np.unravel_index(s, data.shape)
+            hs.Fill( s[0], s[1])
+    if len(doubles) >0 :
+        for d in doubles:
+            d = np.unravel_index(d, data.shape)
+            hd.Fill( d[0], d[1])
+            
+            
+    
+despike = DRSSpikes(user_action = spikecallback)
+
+def mars_spikes( data ):
+    """
+    should search for spikes, just as it is implemented in 
+    mcore/DrsCalib.h in DrsCalib::RemoveSpikes
+    static void RemoveSpikes(float *vec, uint32_t roi)
+    {
+        if (roi<4)
+            return;
+
+        for (size_t ch=0; ch<1440; ch++)
+        {
+            float *p = vec + ch*roi;
+
+            for (size_t i=1; i<roi-2; i++)
+            {
+                if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
+                {
+                    p[i] = (p[i-1]+p[i+1])/2;
+                }
+
+                if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
+                {
+                    p[i] = (p[i-1]+p[i+2])/2;
+                    p[i+1] = p[i];
+                }
+            }
+        }
+    }    
+        
+    """
+
+    #: list, conaining the (chid,slice) tuple of the single spike positions
+    singles = []
+    #: list, conaining the (chid,slice) tuple of the 1st double spike positions
+    doubles = []
+    
+    for chid, pdata in enumerate(data):
+        single_cand = np.where( np.diff(pdata[:-1]) > 25)[0]
+        for cand in single_cand:
+            if -np.diff(pdata[1:])[cand] > 25:
+                singles.append( (chid, cand) )
+        
+        double_cand = np.where( np.diff(pdata[:-1]) > 22 )[0]
+        for cand in double_cand:
+            if cand+1 < len(np.diff(pdata[1:])):
+                if abs(-np.diff(pdata[1:])[cand])<4 and -np.diff(pdata[1:])[cand+1] > 22:
+                    doubles.append( (chid, cand) )
+                
+    
+    return singles, doubles        
+
+
+
+
+for event in run:
+    data = event['data']
+    s, d = mars_spikes(data)
+    sc = event['start_cells']
+    bt = event['board_times']
+    
+    if len(s) >0 :
+        chid.value  = s[0][0]
+        startcell.value = sc[ chid.value ]
+        number_of_singles.value = len(s)
+        for i in s:
+            log = i[1]
+            phys  = (startcell.value+log)%1024
+            
+            hs_log_chid.Fill( chid, log)
+            hs_phys_log.Fill( phys, log)
+            hs_phys_sc.Fill( stace, phys)
+            hs_log_sc.Fill( stace, log)
+            
+            
+    if len(d) >0 :
+        for i in d:
+            log   = i[1]
+            chid  = i[0]
+            stace = sc[ chid ]
+            phys  = (stace+log)%1024
+
+            hd_log_chid.Fill( chid, log)
+            hd_phys_log.Fill( phys, log)
+            hd_phys_sc.Fill( stace, phys)
+            hd_log_sc.Fill( stace, log)
+            hd_log_sc_pp[chid].Fill(stace, log)
+    
+    
+    
+    if event['event_id'] % 100 == 0:
+        print event['event_id']
+        canv.cd(1)
+        hs_log_chid.Draw("COLZ")
+        canv.cd(2)
+        hd_log_chid.Draw("COLZ")
+        canv.cd(3)
+        hs_phys_log.Draw("COLZ")
+        canv.cd(4)
+        hd_phys_log.Draw("COLZ")
+    
+        canv.cd(5)
+        hs_phys_sc.Draw("COLZ")
+        canv.cd(6)
+        hd_phys_sc.Draw("COLZ")
+        canv.cd(7)
+        hs_log_sc.Draw("COLZ")
+        canv.cd(8)
+        hd_log_sc.Draw("COLZ")
+        canv.Modified()
+        canv.Update()
+        
+        
+
+    if event['event_id'] % 1000 == 0:
+        print event['event_id']
+        hs_log_chid.Write('',TObject.kOverwrite)
+        hd_log_chid.Write('',TObject.kOverwrite)
+        hs_phys_log.Write('',TObject.kOverwrite)
+        hd_phys_log.Write('',TObject.kOverwrite)
+        hs_phys_sc.Write('',TObject.kOverwrite)
+        hd_phys_sc.Write('',TObject.kOverwrite)
+        hs_log_sc.Write('',TObject.kOverwrite)
+        hd_log_sc.Write('',TObject.kOverwrite)
+        
+        rootfile.cd('per_chid')
+        for h in hd_log_sc_pp:
+            h.Write('',TObject.kOverwrite)
+        rootfile.cd('..')
+    
+hs_log_chid.Write('',TObject.kOverwrite)
+hd_log_chid.Write('',TObject.kOverwrite)
+hs_phys_log.Write('',TObject.kOverwrite)
+hd_phys_log.Write('',TObject.kOverwrite)
+hs_phys_sc.Write('',TObject.kOverwrite)
+hd_phys_sc.Write('',TObject.kOverwrite)
+hs_log_sc.Write('',TObject.kOverwrite)
+hd_log_sc.Write('',TObject.kOverwrite)
+        
+rootfile.cd('per_chid')
+for h in hd_log_sc_pp:
+    h.Write()
+rootfile.cd('..')    
+rootfile.Close()
+
+
+def book_histos():
+    
+    histograms = {}
+    histograms['overview'] = {}
+    ov = histograms['overview']
+    histograms['per_pixel'] = {}
+    pp = histograms['per_pixel']
+    
+    ov['hs_log_chid'] = TH2F("hs_log_chid",
+        "singles",
+        1440, -0.5, 1439.5, 
+        300, -0.5, 299.5)
+    ov['hs_log_chid'].GetXaxis().SetTitle("cont hardware ID (CHID)")
+    ov['hs_log_chid'].GetYaxis().SetTitle("spike position in the logical pipeline")
+
+    ov['hd_log_chid'] = TH2F("hd_log_chid",
+        "doubles",
+        1440, -0.5, 1439.5, 
+        300, -0.5, 299.5)
+    ov['hd_log_chid'].GetXaxis().SetTitle("cont hardware ID (CHID)")
+    ov['hd_log_chid'].GetYaxis().SetTitle("spike position in the logical pipeline")
+
+
+    ov['hs_phys_log'] = TH2F("hs_phys_log",
+        "singles",
+        1024, -0.5, 1023.5, 
+        300, -0.5, 299.5)
+    ov['hs_phys_log'].GetXaxis().SetTitle("spike position in the physical pipeline")
+    ov['hs_phys_log'].GetYaxis().SetTitle("spike position in the logical pipeline")
+
+    ov['hd_phys_log'] = TH2F("hd_phys_log",
+        "doubles",
+        1024, -0.5, 1023.5, 
+        300, -0.5, 299.5)
+    ov['hd_phys_log'].GetXaxis().SetTitle("spike position in the physical pipeline")
+    ov['hd_phys_log'].GetYaxis().SetTitle("spike position in the logical pipeline")
+
+
+    ov['hs_phys_sc'] = TH2F("hs_phys_sc",
+        "singles",
+        1024, -0.5, 1023.5, 
+        1024, -0.5, 1023.5)
+    ov['hs_phys_sc'].GetXaxis().SetTitle("DRS startcell")
+    ov['hs_phys_sc'].GetYaxis().SetTitle("spike position in the physical pipeline")
+    
+    ov['hd_phys_sc'] = TH2F("hd_phys_sc",
+        "doubles",
+        1024, -0.5, 1023.5, 
+        1024, -0.5, 1023.5)
+    ov['hd_phys_sc'].GetXaxis().SetTitle("DRS startcell")
+    ov['hd_phys_sc'].GetYaxis().SetTitle("spike position in the physical pipeline")
+
+
+
+    ov['hs_log_sc'] = TH2F("hs_log_sc",
+        "singles",
+        1024, -0.5, 1023.5, 
+        1024, -0.5, 1023.5)
+    ov['hs_log_sc'].GetXaxis().SetTitle("DRS startcell")
+    ov['hs_log_sc'].GetYaxis().SetTitle("spike position in the logical pipeline")
+    
+    ov['hd_log_sc'] = TH2F("hd_log_sc",
+        "doubles",
+        1024, -0.5, 1023.5, 
+        1024, -0.5, 1023.5)
+    ov['hd_log_sc'].GetXaxis().SetTitle("DRS startcell")
+    ov['hd_log_sc'].GetYaxis().SetTitle("spike position in the logical pipeline")
+
+
+    hd_log_sc_pp = []
+    for chid in range(1440):
+        hd_log_sc_pp.append(
+            TH2F("chid_"+str(chid),
+                "doubles - log pipe vs startcell - chid:" +str(chid),
+                1024, -0.5, 1023.5, 
+                300, -0.5, 299.5)
+            )
+            
+
+
