source: fact/tools/pyscripts/sandbox/dneise/spikes/ananew.py@ 14788

Last change on this file since 14788 was 14180, checked in by neise, 12 years ago
initial
  • Property svn:executable set to *
File size: 4.3 KB
Line 
1#!/usr/bin/python -tt
2
3from pyfact import RawData
4from drs_spikes import DRSSpikes
5import sys
6import numpy as np
7from ctypes import *
8import os
9import os.path as path
10from ROOT import TFile, TCanvas, TH2F, TTree, TStyle, TObject
11import time
12
13dfn = sys.argv[1]
14cfn = sys.argv[2]
15run = RawData( dfn, cfn)
16
17
18file_base_name = path.splitext(path.basename(dfn))[0]
19root_filename = '/home_nfs/isdc/neise/' + file_base_name + '_' + time.strftime('%Y%m%d_%H%M%S') +'_spikeana.root'
20
21rootfile = TFile(root_filename, "RECREATE")
22
23fbn = path.splitext(file_base_name)[0]
24fbn = fbn[0:8]+fbn[9:12]
25
26
27baum = TTree('spiketree', 'spike ana tree')
28
29# prepare some vars for the tree
30chid = c_int(0)
31startcell = c_int(0)
32number_of_singles = c_int(0)
33position_of_spikes_in_logical_pipeline = np.zeros(100, np.int32)
34position_of_spikes_in_physical_pipeline = np.zeros(100, np.int32)
35time_to_previous = c_int(0)
36event_id = c_int(0)
37file_id = c_long(long(fbn))
38
39baum.Branch('file',file_id,'file/L')
40baum.Branch('event',event_id,'event/I')
41baum.Branch('chid',chid,'chid/I')
42baum.Branch('sc',startcell,'sc/I')
43baum.Branch('n',number_of_singles,'n/I')
44baum.Branch('logpos',position_of_spikes_in_logical_pipeline,'logpos[n]/I')
45baum.Branch('physpos',position_of_spikes_in_physical_pipeline,'physpos[n]/I')
46baum.Branch('time',time_to_previous,'time/I')
47
48
49def spikecallback(candidates, singles, doubles, data, ind):
50 if len(singles) >0 :
51 for s in singles:
52 s = np.unravel_index(s, data.shape)
53 hs.Fill( s[0], s[1])
54 if len(doubles) >0 :
55 for d in doubles:
56 d = np.unravel_index(d, data.shape)
57 hd.Fill( d[0], d[1])
58
59
60
61despike = DRSSpikes(user_action = spikecallback)
62
63def mars_spikes( data ):
64 """
65 should search for spikes, just as it is implemented in
66 mcore/DrsCalib.h in DrsCalib::RemoveSpikes
67 static void RemoveSpikes(float *vec, uint32_t roi)
68 {
69 if (roi<4)
70 return;
71
72 for (size_t ch=0; ch<1440; ch++)
73 {
74 float *p = vec + ch*roi;
75
76 for (size_t i=1; i<roi-2; i++)
77 {
78 if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
79 {
80 p[i] = (p[i-1]+p[i+1])/2;
81 }
82
83 if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
84 {
85 p[i] = (p[i-1]+p[i+2])/2;
86 p[i+1] = p[i];
87 }
88 }
89 }
90 }
91
92 """
93
94 #: list, conaining the (chid,slice) tuple of the single spike positions
95 singles = []
96 #: list, conaining the (chid,slice) tuple of the 1st double spike positions
97 doubles = []
98
99 for chid, pdata in enumerate(data):
100 single_cand = np.where( np.diff(pdata[:-1]) > 25)[0]
101 for cand in single_cand:
102 if -np.diff(pdata[1:])[cand] > 25:
103 singles.append( (chid, cand) )
104
105 double_cand = np.where( np.diff(pdata[:-1]) > 22 )[0]
106 for cand in double_cand:
107 if cand+1 < len(np.diff(pdata[1:])):
108 if abs(-np.diff(pdata[1:])[cand])<4 and -np.diff(pdata[1:])[cand+1] > 22:
109 doubles.append( (chid, cand) )
110
111
112 return singles, doubles
113
114
115
116
117
118event = run.next()
119board_time = event['board_times'].copy()
120
121for event in run:
122 event_id.value = event['event_id']
123
124
125 data = event['data']
126 singles, doubles = mars_spikes(data)
127 sc = event['start_cells']
128 bt_old = board_time.copy()
129 board_time = event['board_times'].copy()
130 # print fbn, event['event_id']
131
132 ss = []
133 for i in range(1440):
134 ss.append([])
135 for s in singles:
136 ss[s[0]].append(s[1])
137
138 for pixel_id,pixel in enumerate(ss):
139 if len(pixel) >0 :
140 chid.value = pixel_id
141 startcell.value = sc[ chid.value ]
142 number_of_singles.value = len(pixel)
143 time_to_previous.value = board_time[pixel_id/36]-bt_old[pixel_id/36]
144 for num,s in enumerate(pixel):
145 log = s
146 phys = (startcell.value+log)%1024
147 position_of_spikes_in_logical_pipeline[num] = log
148 position_of_spikes_in_physical_pipeline[num] = phys
149 baum.Fill()
150
151
152baum.Write()
153
154rootfile.Close()
Note: See TracBrowser for help on using the repository browser.