source: fact/tools/pyscripts/sandbox/dneise/spikes/ana_old.py@ 16858

Last change on this file since 16858 was 14171, checked in by neise, 12 years ago
working
  • Property svn:executable set to *
File size: 8.1 KB
Line 
1#!/usr/bin/python -tti
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
23singlebaum = TTree('singlebaum', 'spike ana tree')
24
25# prepare some vars for the tree
26chid = c_int(0)
27startcell = c_int(0)
28number_of_singles = c_int(0)
29position_of_spikes_in_logical_pipeline = np.zeros(100, np.int32)
30position_of_spikes_in_physical_pipeline = np.zeros(100, np.int32)
31time_to_previous = c_long(0)
32
33singlebaum.Branch('chid',chid,'chid/I')
34singlebaum.Branch('sc',startcell,'sc/I')
35singlebaum.Branch('n',number_of_singles,'n/I')
36singlebaum.Branch('logpos',position_of_spikes_in_logical_pipeline,'logpos[n]/I')
37singlebaum.Branch('physpos',position_of_spikes_in_physical_pipeline,'physpos[n]/I')
38singlebaum.Branch('time',time_to_previous,'time/I')
39
40
41def spikecallback(candidates, singles, doubles, data, ind):
42 if len(singles) >0 :
43 for s in singles:
44 s = np.unravel_index(s, data.shape)
45 hs.Fill( s[0], s[1])
46 if len(doubles) >0 :
47 for d in doubles:
48 d = np.unravel_index(d, data.shape)
49 hd.Fill( d[0], d[1])
50
51
52
53despike = DRSSpikes(user_action = spikecallback)
54
55def mars_spikes( data ):
56 """
57 should search for spikes, just as it is implemented in
58 mcore/DrsCalib.h in DrsCalib::RemoveSpikes
59 static void RemoveSpikes(float *vec, uint32_t roi)
60 {
61 if (roi<4)
62 return;
63
64 for (size_t ch=0; ch<1440; ch++)
65 {
66 float *p = vec + ch*roi;
67
68 for (size_t i=1; i<roi-2; i++)
69 {
70 if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
71 {
72 p[i] = (p[i-1]+p[i+1])/2;
73 }
74
75 if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
76 {
77 p[i] = (p[i-1]+p[i+2])/2;
78 p[i+1] = p[i];
79 }
80 }
81 }
82 }
83
84 """
85
86 #: list, conaining the (chid,slice) tuple of the single spike positions
87 singles = []
88 #: list, conaining the (chid,slice) tuple of the 1st double spike positions
89 doubles = []
90
91 for chid, pdata in enumerate(data):
92 single_cand = np.where( np.diff(pdata[:-1]) > 25)[0]
93 for cand in single_cand:
94 if -np.diff(pdata[1:])[cand] > 25:
95 singles.append( (chid, cand) )
96
97 double_cand = np.where( np.diff(pdata[:-1]) > 22 )[0]
98 for cand in double_cand:
99 if cand+1 < len(np.diff(pdata[1:])):
100 if abs(-np.diff(pdata[1:])[cand])<4 and -np.diff(pdata[1:])[cand+1] > 22:
101 doubles.append( (chid, cand) )
102
103
104 return singles, doubles
105
106
107
108
109for event in run:
110 data = event['data']
111 s, d = mars_spikes(data)
112 sc = event['start_cells']
113 bt = event['board_times']
114
115 if len(s) >0 :
116 chid.value = s[0][0]
117 startcell.value = sc[ chid.value ]
118 number_of_singles.value = len(s)
119 for i in s:
120 log = i[1]
121 phys = (startcell.value+log)%1024
122
123 hs_log_chid.Fill( chid, log)
124 hs_phys_log.Fill( phys, log)
125 hs_phys_sc.Fill( stace, phys)
126 hs_log_sc.Fill( stace, log)
127
128
129 if len(d) >0 :
130 for i in d:
131 log = i[1]
132 chid = i[0]
133 stace = sc[ chid ]
134 phys = (stace+log)%1024
135
136 hd_log_chid.Fill( chid, log)
137 hd_phys_log.Fill( phys, log)
138 hd_phys_sc.Fill( stace, phys)
139 hd_log_sc.Fill( stace, log)
140 hd_log_sc_pp[chid].Fill(stace, log)
141
142
143
144 if event['event_id'] % 100 == 0:
145 print event['event_id']
146 canv.cd(1)
147 hs_log_chid.Draw("COLZ")
148 canv.cd(2)
149 hd_log_chid.Draw("COLZ")
150 canv.cd(3)
151 hs_phys_log.Draw("COLZ")
152 canv.cd(4)
153 hd_phys_log.Draw("COLZ")
154
155 canv.cd(5)
156 hs_phys_sc.Draw("COLZ")
157 canv.cd(6)
158 hd_phys_sc.Draw("COLZ")
159 canv.cd(7)
160 hs_log_sc.Draw("COLZ")
161 canv.cd(8)
162 hd_log_sc.Draw("COLZ")
163 canv.Modified()
164 canv.Update()
165
166
167
168 if event['event_id'] % 1000 == 0:
169 print event['event_id']
170 hs_log_chid.Write('',TObject.kOverwrite)
171 hd_log_chid.Write('',TObject.kOverwrite)
172 hs_phys_log.Write('',TObject.kOverwrite)
173 hd_phys_log.Write('',TObject.kOverwrite)
174 hs_phys_sc.Write('',TObject.kOverwrite)
175 hd_phys_sc.Write('',TObject.kOverwrite)
176 hs_log_sc.Write('',TObject.kOverwrite)
177 hd_log_sc.Write('',TObject.kOverwrite)
178
179 rootfile.cd('per_chid')
180 for h in hd_log_sc_pp:
181 h.Write('',TObject.kOverwrite)
182 rootfile.cd('..')
183
184hs_log_chid.Write('',TObject.kOverwrite)
185hd_log_chid.Write('',TObject.kOverwrite)
186hs_phys_log.Write('',TObject.kOverwrite)
187hd_phys_log.Write('',TObject.kOverwrite)
188hs_phys_sc.Write('',TObject.kOverwrite)
189hd_phys_sc.Write('',TObject.kOverwrite)
190hs_log_sc.Write('',TObject.kOverwrite)
191hd_log_sc.Write('',TObject.kOverwrite)
192
193rootfile.cd('per_chid')
194for h in hd_log_sc_pp:
195 h.Write()
196rootfile.cd('..')
197rootfile.Close()
198
199
200def book_histos():
201
202 histograms = {}
203 histograms['overview'] = {}
204 ov = histograms['overview']
205 histograms['per_pixel'] = {}
206 pp = histograms['per_pixel']
207
208 ov['hs_log_chid'] = TH2F("hs_log_chid",
209 "singles",
210 1440, -0.5, 1439.5,
211 300, -0.5, 299.5)
212 ov['hs_log_chid'].GetXaxis().SetTitle("cont hardware ID (CHID)")
213 ov['hs_log_chid'].GetYaxis().SetTitle("spike position in the logical pipeline")
214
215 ov['hd_log_chid'] = TH2F("hd_log_chid",
216 "doubles",
217 1440, -0.5, 1439.5,
218 300, -0.5, 299.5)
219 ov['hd_log_chid'].GetXaxis().SetTitle("cont hardware ID (CHID)")
220 ov['hd_log_chid'].GetYaxis().SetTitle("spike position in the logical pipeline")
221
222
223 ov['hs_phys_log'] = TH2F("hs_phys_log",
224 "singles",
225 1024, -0.5, 1023.5,
226 300, -0.5, 299.5)
227 ov['hs_phys_log'].GetXaxis().SetTitle("spike position in the physical pipeline")
228 ov['hs_phys_log'].GetYaxis().SetTitle("spike position in the logical pipeline")
229
230 ov['hd_phys_log'] = TH2F("hd_phys_log",
231 "doubles",
232 1024, -0.5, 1023.5,
233 300, -0.5, 299.5)
234 ov['hd_phys_log'].GetXaxis().SetTitle("spike position in the physical pipeline")
235 ov['hd_phys_log'].GetYaxis().SetTitle("spike position in the logical pipeline")
236
237
238 ov['hs_phys_sc'] = TH2F("hs_phys_sc",
239 "singles",
240 1024, -0.5, 1023.5,
241 1024, -0.5, 1023.5)
242 ov['hs_phys_sc'].GetXaxis().SetTitle("DRS startcell")
243 ov['hs_phys_sc'].GetYaxis().SetTitle("spike position in the physical pipeline")
244
245 ov['hd_phys_sc'] = TH2F("hd_phys_sc",
246 "doubles",
247 1024, -0.5, 1023.5,
248 1024, -0.5, 1023.5)
249 ov['hd_phys_sc'].GetXaxis().SetTitle("DRS startcell")
250 ov['hd_phys_sc'].GetYaxis().SetTitle("spike position in the physical pipeline")
251
252
253
254 ov['hs_log_sc'] = TH2F("hs_log_sc",
255 "singles",
256 1024, -0.5, 1023.5,
257 1024, -0.5, 1023.5)
258 ov['hs_log_sc'].GetXaxis().SetTitle("DRS startcell")
259 ov['hs_log_sc'].GetYaxis().SetTitle("spike position in the logical pipeline")
260
261 ov['hd_log_sc'] = TH2F("hd_log_sc",
262 "doubles",
263 1024, -0.5, 1023.5,
264 1024, -0.5, 1023.5)
265 ov['hd_log_sc'].GetXaxis().SetTitle("DRS startcell")
266 ov['hd_log_sc'].GetYaxis().SetTitle("spike position in the logical pipeline")
267
268
269 hd_log_sc_pp = []
270 for chid in range(1440):
271 hd_log_sc_pp.append(
272 TH2F("chid_"+str(chid),
273 "doubles - log pipe vs startcell - chid:" +str(chid),
274 1024, -0.5, 1023.5,
275 300, -0.5, 299.5)
276 )
277
278
279
Note: See TracBrowser for help on using the repository browser.