source: fact/tools/pyscripts/new_pyfact/pyfact.py@ 17972

Last change on this file since 17972 was 17972, checked in by dneise, 10 years ago
added convenience class AuxFile
  • Property svn:executable set to *
File size: 12.9 KB
Line 
1#!/usr/bin/python -tt
2#
3# Werner Lustermann, Dominik Neise
4# ETH Zurich, TU Dortmund
5#
6import os
7import sys
8import numpy as np
9import ROOT
10import pylab
11
12from scipy.signal import firwin, iirdesign, lfilter
13
14########## BUILDING OF THE SHARED OBJECT FILES ###############################
15if __name__ == '__main__' and len(sys.argv) > 1 and 'build' in sys.argv[1]:
16 ROOT.gSystem.AddLinkedLibs("-lz")
17 root_make_string = ROOT.gSystem.GetMakeSharedLib()
18 if not "-std=c++0x" in root_make_string:
19 make_string = root_make_string.replace('$Opt', '$Opt -std=c++0x -D HAVE_ZLIB')
20 ROOT.gSystem.SetMakeSharedLib(make_string)
21 ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/izstream.h+O")
22 ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/fits.h+O")
23 ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/zfits.h+O")
24 ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/factfits.h+O")
25 if not "-std=c++0x" in root_make_string:
26 make_string = root_make_string.replace('$Opt', "$Opt -std=c++0x -D HAVE_ZLIB -D'PACKAGE_NAME=\"PACKAGE_NAME\"' "
27 "-D'PACKAGE_VERSION=\"PACKAGE_VERSION\"' -D'REVISION=\"REVISION\"' ")
28 ROOT.gSystem.SetMakeSharedLib(make_string)
29 ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/DrsCalib.h+O")
30
31 ROOT.gInterpreter.GenerateDictionary("map<string,fits::Entry>","map;string;extern_Mars_mcore/fits.h")
32 ROOT.gInterpreter.GenerateDictionary("pair<string,fits::Entry>","map;string;extern_Mars_mcore/fits.h")
33 ROOT.gInterpreter.GenerateDictionary("map<string,fits::Table::Column>","map;string;extern_Mars_mcore/fits.h")
34 ROOT.gInterpreter.GenerateDictionary("pair<string,fits::Table::Column>","map;string;extern_Mars_mcore/fits.h")
35 ROOT.gInterpreter.GenerateDictionary("vector<DrsCalibrate::Step>","vector;extern_Mars_mcore/DrsCalib.h")
36
37########## USAGE #############################################################
38if __name__ == '__main__' and len(sys.argv) < 3:
39 print """ Usage:
40 ----------------------------------------------------------------------
41 To just build the shared object libs call:
42 python pyfact.py build
43
44 To build (if necessary) and open an example file
45 python -i pyfact.py /path/to/data_file.zfits /path/to/calib_file.drs.fits.gz
46
47 Any of the 3 file 'types': fits.gz zfits fits
48 should be supported.
49
50 To clean all of automatically built files, do something like:
51 rm *.so *.d AutoDict_* extern_Mars_mcore/*.so extern_Mars_mcore/*.d
52 ----------------------------------------------------------------------
53 """
54 sys.exit(1)
55
56######### START OF PYFACT MODULE ############################################
57path = os.path.dirname(os.path.realpath(__file__))
58ROOT.gSystem.Load(path+'/AutoDict_map_string_fits__Entry__cxx.so')
59ROOT.gSystem.Load(path+'/AutoDict_map_string_fits__Table__Column__cxx.so')
60ROOT.gSystem.Load(path+'/AutoDict_pair_string_fits__Entry__cxx.so')
61ROOT.gSystem.Load(path+'/AutoDict_pair_string_fits__Table__Column__cxx.so')
62ROOT.gSystem.Load(path+'/AutoDict_vector_DrsCalibrate__Step__cxx.so')
63ROOT.gSystem.Load(path+'/extern_Mars_mcore/fits_h.so')
64ROOT.gSystem.Load(path+'/extern_Mars_mcore/izstream_h.so')
65ROOT.gSystem.Load(path+'/extern_Mars_mcore/zfits_h.so')
66ROOT.gSystem.Load(path+'/extern_Mars_mcore/factfits_h.so')
67ROOT.gSystem.Load(path+'/extern_Mars_mcore/DrsCalib_h.so')
68del path
69
70class Fits( object ):
71 """ General FITS file access
72
73 Wrapper for factfits class from Mars/mcore/factfits.h
74 (factfits might not be suited to read any FITS file out there, but certainly
75 all FACT FITS files can be read using this class)
76 """
77 __module__ = 'pyfact'
78 def __init__(self, path):
79 """
80 """
81 if not os.path.exists(path):
82 raise IOError(path+' was not found')
83 self.f = ROOT.factfits(path)
84 self._make_header()
85 self._setup_columns()
86
87 def _make_header(self):
88 """
89 """
90 str_to_bool = { 'T':True, 'F':False}
91 type_conversion = { 'I' : int, 'F' : float, 'T' : str, 'B' : str_to_bool.__getitem__}
92
93 self.header = {}
94 self.header_comments = {}
95 for key,entry in self.f.GetKeys():
96 try:
97 self.header[key] = type_conversion[entry.type](entry.value)
98 except KeyError:
99 raise IOError("Error: entry type unknown.\n Is %s, but should be one of: [I,F,T,B]" % (entry.type) )
100 self.header_comments[key] = entry.comment
101
102 def _setup_columns(self):
103 """
104 """
105 col_type_to_np_type_map = { 'L' : 'b1', 'A' : 'a1', 'B' : 'i1',
106 'I' : 'i2', 'J' : 'i4', 'K' : 'i8', 'E' : 'f4', 'D' : 'f8'}
107 self.cols = {}
108 for key,col in self.f.GetColumns():
109 self.cols[key] = np.zeros(col.num, col_type_to_np_type_map[col.type])
110 if col.num != 0:
111 self.f.SetPtrAddress(key, self.cols[key])
112
113 def __iter__(self):
114 return self
115
116 def next(self, row=None):
117 """
118 """
119 if row is None:
120 if self.f.GetNextRow() == False:
121 raise StopIteration
122 else:
123 row = int(row)
124 if self.f.GetRow(row) == False:
125 raise StopIteration
126 return self
127
128class AuxFile( Fits ):
129 """ easy(?) access to FACT aux files
130 """
131 __module__ = 'pyfact'
132 def __init__(self, path, verbose=False):
133 self._verbose = verbose
134 super(AuxFile, self).__init__(path)
135
136 def _setup_columns(self):
137 col_type_to_np_type_map = { 'L' : 'b1', 'A' : 'a1', 'B' : 'i1',
138 'I' : 'i2', 'J' : 'i4', 'K' : 'i8', 'E' : 'f4', 'D' : 'f8'}
139 self._cols = {}
140 self.cols = {}
141 N = self.header['NAXIS2']
142 for key,col in self.f.GetColumns():
143 self._cols[key] = np.zeros(col.num, col_type_to_np_type_map[col.type])
144 self.cols[key] = np.zeros((N, col.num), col_type_to_np_type_map[col.type])
145 if col.num != 0:
146 self.f.SetPtrAddress(key, self._cols[key])
147
148
149 for i,row in enumerate(self):
150 if self._verbose:
151 try:
152 step = int(self._verbose)
153 except:
154 step = 10
155 if i % step == 0:
156 print "reading line", i
157
158 for key in self._cols:
159 self.cols[key][i,:] = self._cols[key]
160
161
162
163
164class RawData( Fits ):
165 """ Special raw data FITS file access (with DRS4 calibration)
166
167 During iteration the C++ method DrsCalibration::Apply is being called.
168 """
169 __module__='pyfact'
170 def __init__(self, data_path, calib_path):
171 """ -constructor-
172 *data_path* : fits or fits.gz file of the data including the path
173 *calib_path* : fits or fits.gz file containing DRS calibration data
174 """
175 super(RawData, self).__init__(data_path)
176 self.cols['CalibData'] = np.zeros( self.cols['Data'].shape, np.float32)
177 self.cols['CalibData2D'] = self.cols['CalibData'].reshape( self.header['NPIX'], -1)
178 if not self.cols['CalibData2D'].base is self.cols['CalibData']:
179 print "Error seomthing went wrong!"
180 self.drs_calibration = ROOT.DrsCalibration()
181 self.drs_calibration.ReadFitsImp( calib_path )
182
183 self.drs_calibrate = ROOT.DrsCalibrate()
184 self.list_of_previous_start_cells = []
185
186 def __iter__(self):
187 """ iterator """
188 return self
189
190 def next(self, row=None):
191 """
192 """
193 super(RawData, self).next(row)
194
195 self.drs_calibration.Apply( self.cols['CalibData'],
196 self.cols['Data'],
197 self.cols['StartCellData'],
198 self.header['NROI'])
199
200 for previous_start_cells in self.list_of_previous_start_cells:
201 self.drs_calibrate.CorrectStep(
202 self.cols['CalibData'],
203 self.header['NPIX'],
204 self.header['NROI'],
205 previous_start_cells,
206 self.cols['StartCellData'],
207 self.header['NROI']+10)
208 self.drs_calibrate.CorrectStep(
209 self.cols['CalibData'],
210 self.header['NPIX'],
211 self.header['NROI'],
212 previous_start_cells,
213 self.cols['StartCellData'],
214 3)
215 self.list_of_previous_start_cells.append(self.cols['StartCellData'])
216 if len(self.list_of_previous_start_cells) > 5:
217 self.list_of_previous_start_cells.pop(0)
218
219 for ch in range(self.header['NPIX']):
220 self.drs_calibrate.RemoveSpikes3(self.cols['CalibData2D'][ch], self.header['NROI'])
221
222 return self
223
224
225
226if __name__ == '__main__':
227 """ Example """
228 print "Example for calibrated raw-file"
229 f = RawData(sys.argv[1], sys.argv[2])
230
231 print "number of events:", f.header['NAXIS2']
232 print "date of observation:", f.header['DATE']
233 print "The files has these cols:", f.cols.keys()
234
235 for counter,row in enumerate(f):
236 print "Event Id:", row.cols['EventNum']
237 print "shape of column 'StartCellData'", row.cols['StartCellData'].shape
238 print "dtype of column 'Data'", row.cols['StartCellData'].dtype
239 if counter > 3:
240 break
241 # get next row
242 f.next()
243 print "Event Id:", f.cols['EventNum']
244 # get another row
245 f.next(10)
246 print "Event Id:", f.cols['EventNum']
247 # Go back again
248 f.next(3)
249 print "Event Id:", f.cols['EventNum']
250
251 import matplotlib.pyplot as plt
252 plt.ion()
253 """
254 for i in range(f.header['NPIX']):
255 plt.cla()
256 plt.title("Event %d"%(f.cols['EventNum']))
257 plt.plot(f.cols['CalibData2D'][i], '.:', label='pixel %d'%i)
258 plt.legend()
259 answer = raw_input('anykey for next pixel; "q" to quit. :')
260 if 'q' in answer.lower():
261 break
262 """
263
264
265
266 from scipy import signal
267 from scipy.signal import freqs, iirfilter
268 #b, a = iirfilter(4, [1, 10], 1, 60, analog=True, ftype='cheby1')
269 #b, a = signal.iirfilter(17, [50, 200], rs=60, btype='band',
270 # analog=True, ftype='cheby2')
271 #w, h = signal.freqs(b, a, 1000)
272 #w, h = freqs(b, a, worN=np.logspace(-1, 2, 1000))
273
274 b, a = iirfilter(N=4, # the order of the filter
275 Wn = [0.15, 0.25], # critical frequencies
276 rp = 1, # maximum ripple in passband
277 rs = 60, # minimum attenuation in stopband
278 btype='bandstop', # {'bandpass', 'lowpass', 'highpass', 'bandstop'}
279 #analog=True,
280 ftype='cheby1')
281 #w, h = freqs(b, a, worN=np.logspace(-1, 5, 1000))
282 w, h = freqs(b, a)
283 plt.figure()
284 plt.semilogx(w, 20 * np.log10(abs(h)))
285 plt.xlabel('Frequency')
286 plt.ylabel('Amplitude response [dB]')
287 plt.grid()
288 plt.show()
289
290
291
292
293
294 fig1, (ax11,ax12) = plt.subplots(nrows=2)
295 fig2, (ax21,ax22) = plt.subplots(nrows=2)
296
297 flt = firwin(13, [0.15,0.25])
298 t = np.arange(f.cols['CalibData2D'].shape[1])/2.*1e-9
299 fig3, (ax31, ax32) = plt.subplots(nrows=2)
300 ax31.plot(flt)
301 #w, h = freqs(flt, [1.0])
302 w, h = freqs(flt, [1.0], worN=np.logspace(-1, 5, 1000))
303 ax32.semilogx(w, 20 * np.log10(abs(h)))
304 ax32.set_xlabel('Frequency')
305 ax32.set_ylabel('Amplitude response [dB]')
306 ax32.grid()
307
308
309
310 for i in range(f.header['NPIX']):
311 data = f.cols['CalibData2D'][i]
312 fil_data = lfilter(flt, [1.0], data)
313
314 ax11.cla()
315 ax12.cla()
316 ax11.set_title("Event %d"%(f.cols['EventNum']))
317 ax11.plot(t, data, '.:', label='pixel %d'%i)
318 ax11.legend()
319 ax11.set_ylim(-20,40)
320 ax12.plot(t[:fil_data.size], fil_data, '.:', label='pixel %d'%i)
321 ax12.legend()
322 ax12.set_ylim(-20,40)
323
324 """
325 Pxx, freqs, t, plot = pylab.specgram( f.cols['CalibData2D'][i],
326 NFFT=32, Fs=2000000000,
327 detrend=pylab.detrend_none,
328 window=pylab.window_hanning,
329 noverlap=16)
330 """
331 ax21.cla()
332 ax22.cla()
333 dat, freqs, bins, im = ax21.specgram(data,
334 NFFT=64, Fs=2000000000,
335 detrend=pylab.detrend_none,
336 window=pylab.window_hanning,
337 noverlap=60,
338 interpolation='nearest')
339 ax21.axis('tight')
340 ax21.set_ylabel("frequency [Hz]")
341 ax21.set_xlabel("time [s]")
342
343 dat2, freqs2, bins2, im2 = ax22.specgram(data,
344 NFFT=64, Fs=2000000000,
345 detrend=pylab.detrend_none,
346 window=pylab.window_hanning,
347 noverlap=60,
348 interpolation='nearest')
349 ax22.axis('tight')
350 ax22.set_ylabel("frequency [Hz]")
351 ax22.set_xlabel("time [s]")
352
353 plt.show()
354
355 answer = raw_input('anykey for next pixel; "q" to quit. :')
356 if 'q' in answer.lower():
357 break
Note: See TracBrowser for help on using the repository browser.