source: fact/tools/pyscripts/sandbox/dneise/spikes/ana.py

Last change on this file was 14175, checked in by neise, 12 years ago
can now create trees
  • Property svn:executable set to *
File size: 4.0 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
23baum = TTree('spiketree', '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
33
34baum.Branch('chid',chid,'chid/I')
35baum.Branch('sc',startcell,'sc/I')
36baum.Branch('n',number_of_singles,'n/I')
37baum.Branch('logpos',position_of_spikes_in_logical_pipeline,'logpos[n]/I')
38baum.Branch('physpos',position_of_spikes_in_physical_pipeline,'physpos[n]/I')
39baum.Branch('time',time_to_previous,'time/I')
40
41
42def spikecallback(candidates, singles, doubles, data, ind):
43 if len(singles) >0 :
44 for s in singles:
45 s = np.unravel_index(s, data.shape)
46 hs.Fill( s[0], s[1])
47 if len(doubles) >0 :
48 for d in doubles:
49 d = np.unravel_index(d, data.shape)
50 hd.Fill( d[0], d[1])
51
52
53
54despike = DRSSpikes(user_action = spikecallback)
55
56def mars_spikes( data ):
57 """
58 should search for spikes, just as it is implemented in
59 mcore/DrsCalib.h in DrsCalib::RemoveSpikes
60 static void RemoveSpikes(float *vec, uint32_t roi)
61 {
62 if (roi<4)
63 return;
64
65 for (size_t ch=0; ch<1440; ch++)
66 {
67 float *p = vec + ch*roi;
68
69 for (size_t i=1; i<roi-2; i++)
70 {
71 if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
72 {
73 p[i] = (p[i-1]+p[i+1])/2;
74 }
75
76 if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
77 {
78 p[i] = (p[i-1]+p[i+2])/2;
79 p[i+1] = p[i];
80 }
81 }
82 }
83 }
84
85 """
86
87 #: list, conaining the (chid,slice) tuple of the single spike positions
88 singles = []
89 #: list, conaining the (chid,slice) tuple of the 1st double spike positions
90 doubles = []
91
92 for chid, pdata in enumerate(data):
93 single_cand = np.where( np.diff(pdata[:-1]) > 25)[0]
94 for cand in single_cand:
95 if -np.diff(pdata[1:])[cand] > 25:
96 singles.append( (chid, cand) )
97
98 double_cand = np.where( np.diff(pdata[:-1]) > 22 )[0]
99 for cand in double_cand:
100 if cand+1 < len(np.diff(pdata[1:])):
101 if abs(-np.diff(pdata[1:])[cand])<4 and -np.diff(pdata[1:])[cand+1] > 22:
102 doubles.append( (chid, cand) )
103
104
105 return singles, doubles
106
107
108
109
110
111event = run.next()
112bt_old = event['board_times'].copy()
113
114for event in run:
115 data = event['data']
116 singles, doubles = mars_spikes(data)
117 sc = event['start_cells']
118 board_time = event['board_times'].copy()
119
120 ss = []
121 for i in range(1440):
122 ss.append([])
123 for s in singles:
124 ss[s[0]].append(s[1])
125
126 for pixel_id,pixel in enumerate(ss):
127 if len(pixel) >0 :
128 chid.value = pixel_id
129 startcell.value = sc[ chid.value ]
130 number_of_singles.value = len(pixel)
131 time.value = board_time[pixel_id/36]-bt_old[pixel_id/36]
132 for num,s in enumerate(pixel):
133 log = s
134 phys = (startcell.value+log)%1024
135 position_of_spikes_in_logical_pipeline[num] = log
136 position_of_spikes_in_physical_pipeline[num] = phys
137 baum.Fill()
138
139 bt_old = event['board_times'].copy()
140
141baum.Write()
142
143rootfile.Close()
Note: See TracBrowser for help on using the repository browser.