1 | #!/usr/bin/python -tt
|
---|
2 | #
|
---|
3 | # Werner Lustermann, Dominik Neise
|
---|
4 | # ETH Zurich, TU Dortmund
|
---|
5 | #
|
---|
6 | from ctypes import *
|
---|
7 | import numpy as np
|
---|
8 | import pprint # for SlowData
|
---|
9 | from scipy import signal
|
---|
10 |
|
---|
11 | # get the ROOT stuff + my shared libs
|
---|
12 | from ROOT import gSystem
|
---|
13 | # factfits_h.so is made from factfits.h and is used to access the data
|
---|
14 | # make sure the location of factfits_h.so is in LD_LIBRARY_PATH.
|
---|
15 | # having it in PYTHONPATH is *not* sufficient
|
---|
16 | gSystem.Load('factfits_h.so')
|
---|
17 | gSystem.Load('calfactfits_h.so')
|
---|
18 | from ROOT import *
|
---|
19 |
|
---|
20 | class RawDataFeeder( object ):
|
---|
21 | """ Wrapper class for RawData class
|
---|
22 | capable of iterating over multiple RawData Files
|
---|
23 | """
|
---|
24 |
|
---|
25 | def __init__(self, filelist):
|
---|
26 | """ *filelist* list of files to iterate over
|
---|
27 | the list should contain tuples, or sublists of two filenames
|
---|
28 | the first should be a data file (\*.fits.gz)
|
---|
29 | the second should be an amplitude calibration file(\*.drs.fits.gz)
|
---|
30 | """
|
---|
31 |
|
---|
32 | self.__module__ = 'pyfact'
|
---|
33 |
|
---|
34 | # sanity check for input
|
---|
35 | if type(filelist) != type(list()):
|
---|
36 | raise TypeError('filelist should be a list')
|
---|
37 | for entry in filelist:
|
---|
38 | if len(entry) != 2:
|
---|
39 | raise TypeError('the entries of filelist should have length == 2')
|
---|
40 | for path in entry:
|
---|
41 | if type(path) != type(str()):
|
---|
42 | raise TypeError('the entries of filelist should be path, i.e. of type str()')
|
---|
43 | #todo check if 'path' is a valid path
|
---|
44 | # else: throw an Exception, or Warning?
|
---|
45 |
|
---|
46 | self.filelist = filelist
|
---|
47 | self._current_RawData = RawData(filelist[0][0], filelist[0][1], return_dict=True)
|
---|
48 | del filelist[0]
|
---|
49 |
|
---|
50 | def __iter__(self):
|
---|
51 | return self
|
---|
52 |
|
---|
53 | def next():
|
---|
54 | """ Method being called by the iterator.
|
---|
55 | Since the RawData Objects are simply looped over, the event_id from the
|
---|
56 | RawData object will not be unique.
|
---|
57 | Each RawData obejct will start with event_id = 1 as usual.
|
---|
58 | """
|
---|
59 | try:
|
---|
60 | return self._current_RawData.next()
|
---|
61 | except StopIteration:
|
---|
62 | # current_RawData was completely processed
|
---|
63 | # delete it (I hope this calls the destructor of the fits file and/or closes it)
|
---|
64 | del self._current_RawData
|
---|
65 | # and remake it, if possible
|
---|
66 | if len(self.filelist) > 0:
|
---|
67 | self._current_RawData = RawData(filelist[0][0], filelist[0][1], return_dict=True)
|
---|
68 | del filelist[0]
|
---|
69 | else:
|
---|
70 | raise
|
---|
71 |
|
---|
72 |
|
---|
73 |
|
---|
74 | class RawData( object ):
|
---|
75 | """ raw data access and calibration
|
---|
76 |
|
---|
77 | class is **iterable**
|
---|
78 |
|
---|
79 | - open raw data file and drs calibration file
|
---|
80 | - performs amplitude calibration
|
---|
81 | - performs baseline substraction if wanted
|
---|
82 | - provides all data in an array:
|
---|
83 | row = number of pixel
|
---|
84 | col = length of region of interest
|
---|
85 |
|
---|
86 | """
|
---|
87 |
|
---|
88 |
|
---|
89 | def __init__(self, data_file_name, calib_file_name,
|
---|
90 | baseline_file_name='',
|
---|
91 | return_dict = True,
|
---|
92 | use_CalFactFits = True,
|
---|
93 | do_calibration = True,
|
---|
94 | user_action_calib=lambda acal_data, data, blm, tom, gm, scells, nroi: None):
|
---|
95 | """ -constructor-
|
---|
96 |
|
---|
97 | - open data file and calibration data file
|
---|
98 | - get basic information about the data in data_file_name
|
---|
99 | - allocate buffers for data access
|
---|
100 |
|
---|
101 | *data_file_name* : fits or fits.gz file of the data including the path
|
---|
102 |
|
---|
103 | *calib_file_name* : fits or fits.gz file containing DRS calibration data
|
---|
104 |
|
---|
105 | *baseline_file_name* : npy file containing the baseline values
|
---|
106 |
|
---|
107 | *return_dict* : this option will be removed in future releases.
|
---|
108 | formerly the next() method returned only a subset of (important) event information,
|
---|
109 | and it was not transparent how to retrieve the other (less important) information.
|
---|
110 | Nowadays next() returns self.__dict__ which contains everything we were able to find in the fits file.
|
---|
111 |
|
---|
112 | *use_CalFactFits* : formerly the DRS amplitude calibration was
|
---|
113 | implemented in python. But for performance reasons this was now moved into
|
---|
114 | a C++ class called CalFactFits. For test purposes, this option can be set to
|
---|
115 | False, but this is not really maintained anymore. If DRS the DRS calibration algorithm is
|
---|
116 | beeing updated in C++ it may not be updated in the python implementation.
|
---|
117 |
|
---|
118 | *do_calibration* : In case *use_CalFactFits* is False, one may choose
|
---|
119 | not to calibrate the data at all, thus safe quite some time.
|
---|
120 | This is imho only needed in case one is interesting in learning something about the
|
---|
121 | calibration algorithm itself.
|
---|
122 |
|
---|
123 | *user_action_calib* : callback function, intended for tests of the DRS calibration algorithm.
|
---|
124 | but since this is not done in the Python regime anymore, this function is never called.
|
---|
125 | (depending on *use_CalFactFits* of course)
|
---|
126 | """
|
---|
127 | self.__module__='pyfact'
|
---|
128 | # manual implementation of default value, but I need to find out
|
---|
129 | # if the user of this class is aware of the new option
|
---|
130 | if return_dict == False:
|
---|
131 | print 'DEPRECATION WARNING:'
|
---|
132 | print 'you are using RawData in a way, which is nor supported anymore.'
|
---|
133 | print ' Please set: return_dict = True, in the __init__ call'
|
---|
134 | self.return_dict = return_dict
|
---|
135 | self.use_CalFactFits = use_CalFactFits
|
---|
136 |
|
---|
137 | self.do_calibration = do_calibration
|
---|
138 |
|
---|
139 | self.data_file_name = data_file_name
|
---|
140 | self.calib_file_name = calib_file_name
|
---|
141 | self.baseline_file_name = baseline_file_name
|
---|
142 |
|
---|
143 | self.user_action_calib = user_action_calib
|
---|
144 |
|
---|
145 | # baseline correction: True / False
|
---|
146 | if len(baseline_file_name) == 0:
|
---|
147 | self.correct_baseline = False
|
---|
148 | else:
|
---|
149 | self.correct_baseline = True
|
---|
150 |
|
---|
151 |
|
---|
152 | # access data file
|
---|
153 | if use_CalFactFits:
|
---|
154 | try:
|
---|
155 | data_file = CalFactFits(data_file_name, calib_file_name)
|
---|
156 | except IOError:
|
---|
157 | print 'problem accessing data file: ', data_file_name
|
---|
158 | raise # stop ! no data
|
---|
159 |
|
---|
160 | #: either CalFactFits object or FactFits object, depending on *use_CalFactFits*
|
---|
161 | self.data_file = data_file
|
---|
162 | #: 1440x300 nparray containing the event data. pixel sorted according to CHID
|
---|
163 | self.data = np.empty( data_file.npix * data_file.nroi, np.float64)
|
---|
164 | data_file.SetNpcaldataPtr(self.data)
|
---|
165 | self.data = self.data.reshape( data_file.npix, data_file.nroi )
|
---|
166 | #: copy of data. here for historical reasons
|
---|
167 | self.acal_data = self.data
|
---|
168 | #: region of interest. (number of DRS slices read).
|
---|
169 | # for FACT data mostly 300. for special runs sometimes 1024.
|
---|
170 | self.nroi = data_file.nroi
|
---|
171 | #: number of Pixel in FACT. should be 1440
|
---|
172 | self.npix = data_file.npix
|
---|
173 | #: the total number of events in the data_file
|
---|
174 | self.nevents = data_file.nevents
|
---|
175 |
|
---|
176 | # Data per event
|
---|
177 | #: starting at 1
|
---|
178 | self.event_id = None
|
---|
179 | #: data=4 ; the rest I don't know by heart .. should be documented here :-)
|
---|
180 | self.trigger_type = None
|
---|
181 | #self.start_cells = None
|
---|
182 | #self.board_times = None
|
---|
183 | #: slice where drs readout started for all DRS chips (160) .. but enlarged to the size of 1440 pixel. thus there are always 9 equal numbers inside.
|
---|
184 | self.start_cells = np.zeros( self.npix, np.int16 )
|
---|
185 | #: each FAD has an onboard clock running from startup time. Currently I don't know the time unit. However this is an array of 40 times, since we have 40 boards.
|
---|
186 | self.board_times = np.zeros( 40, np.int32 )
|
---|
187 |
|
---|
188 | # data_file is a CalFactFits object
|
---|
189 | # data_file.datafile is one of the two FactFits objects hold by a CalFactFits.
|
---|
190 | # sorry for the strange naming ..
|
---|
191 | data_file.datafile.SetPtrAddress('StartCellData', self.start_cells)
|
---|
192 | data_file.datafile.SetPtrAddress('BoardTime', self.board_times)
|
---|
193 |
|
---|
194 | else:
|
---|
195 | try:
|
---|
196 | data_file = FactFits(self.data_file_name)
|
---|
197 | except IOError:
|
---|
198 | print 'problem accessing data file: ', data_file_name
|
---|
199 | raise # stop ! no data
|
---|
200 |
|
---|
201 | self.data_file = data_file
|
---|
202 |
|
---|
203 | # get basic information about the data file
|
---|
204 | self.nroi = data_file.GetUInt('NROI')
|
---|
205 | self.npix = data_file.GetUInt('NPIX')
|
---|
206 | self.nevents = data_file.GetNumRows()
|
---|
207 |
|
---|
208 | # allocate the data memories
|
---|
209 | self.event_id = c_ulong()
|
---|
210 | self.trigger_type = c_ushort()
|
---|
211 | self.data = np.zeros( self.npix * self.nroi, np.int16 ).reshape(self.npix ,self.nroi)
|
---|
212 | self.start_cells = np.zeros( self.npix, np.int16 )
|
---|
213 | self.board_times = np.zeros( 40, np.int32 )
|
---|
214 |
|
---|
215 | # set the pointers to the data++
|
---|
216 | data_file.SetPtrAddress('EventNum', self.event_id)
|
---|
217 | data_file.SetPtrAddress('TriggerType', self.trigger_type)
|
---|
218 | data_file.SetPtrAddress('StartCellData', self.start_cells)
|
---|
219 | data_file.SetPtrAddress('Data', self.data)
|
---|
220 | data_file.SetPtrAddress('BoardTime', self.board_times)
|
---|
221 |
|
---|
222 | # open the calibration file
|
---|
223 | try:
|
---|
224 | calib_file = FactFits(self.calib_file_name)
|
---|
225 | except IOError:
|
---|
226 | print 'problem accessing calibration file: ', calib_file_name
|
---|
227 | raise
|
---|
228 | #: drs calibration file
|
---|
229 | self.calib_file = calib_file
|
---|
230 |
|
---|
231 | baseline_mean = calib_file.GetN('BaselineMean')
|
---|
232 | gain_mean = calib_file.GetN('GainMean')
|
---|
233 | trigger_offset_mean = calib_file.GetN('TriggerOffsetMean')
|
---|
234 |
|
---|
235 | self.Nblm = baseline_mean / self.npix
|
---|
236 | self.Ngm = gain_mean / self.npix
|
---|
237 | self.Ntom = trigger_offset_mean / self.npix
|
---|
238 |
|
---|
239 | self.blm = np.zeros(baseline_mean, np.float32).reshape(self.npix , self.Nblm)
|
---|
240 | self.gm = np.zeros(gain_mean, np.float32).reshape(self.npix , self.Ngm)
|
---|
241 | self.tom = np.zeros(trigger_offset_mean, np.float32).reshape(self.npix , self.Ntom)
|
---|
242 |
|
---|
243 | calib_file.SetPtrAddress('BaselineMean', self.blm)
|
---|
244 | calib_file.SetPtrAddress('GainMean', self.gm)
|
---|
245 | calib_file.SetPtrAddress('TriggerOffsetMean', self.tom)
|
---|
246 | calib_file.GetRow(0)
|
---|
247 |
|
---|
248 | # make calibration constants double, so we never need to roll
|
---|
249 | self.blm = np.hstack((self.blm, self.blm))
|
---|
250 | self.gm = np.hstack((self.gm, self.gm))
|
---|
251 | self.tom = np.hstack((self.tom, self.tom))
|
---|
252 |
|
---|
253 | self.v_bsl = np.zeros(self.npix) # array of baseline values (all ZERO)
|
---|
254 |
|
---|
255 | def __iter__(self):
|
---|
256 | """ iterator """
|
---|
257 | return self
|
---|
258 |
|
---|
259 | def next(self):
|
---|
260 | """ used by __iter__
|
---|
261 |
|
---|
262 | returns self.__dict__
|
---|
263 | """
|
---|
264 | if self.use_CalFactFits:
|
---|
265 | if self.data_file.GetCalEvent() == False:
|
---|
266 | raise StopIteration
|
---|
267 | else:
|
---|
268 | self.event_id = self.data_file.event_id
|
---|
269 | self.trigger_type = self.data_file.event_triggertype
|
---|
270 | #self.start_cells = self.data_file.event_offset
|
---|
271 | #self.board_times = self.data_file.event_boardtimes
|
---|
272 | #self.acal_data = self.data.copy().reshape(self.data_file.npix, self.data_file.nroi)
|
---|
273 | else:
|
---|
274 | if self.data_file.GetNextRow() == False:
|
---|
275 | raise StopIteration
|
---|
276 | else:
|
---|
277 | if self.do_calibration == True:
|
---|
278 | self.calibrate_drs_amplitude()
|
---|
279 |
|
---|
280 | #print 'nevents = ', self.nevents, 'event_id = ', self.event_id.value
|
---|
281 | if self.return_dict:
|
---|
282 | return self.__dict__
|
---|
283 | else:
|
---|
284 | return self.acal_data, self.start_cells, self.trigger_type.value
|
---|
285 |
|
---|
286 | def next_event(self):
|
---|
287 | """ ---- DEPRICATED ----
|
---|
288 |
|
---|
289 | load the next event from disk and calibrate it
|
---|
290 | """
|
---|
291 | if self.use_CalFactFits:
|
---|
292 | self.data_file.GetCalEvent()
|
---|
293 | else:
|
---|
294 | self.data_file.GetNextRow()
|
---|
295 | self.calibrate_drs_amplitude()
|
---|
296 |
|
---|
297 | def calibrate_drs_amplitude(self):
|
---|
298 | """ --- DEPRICATED ---
|
---|
299 |
|
---|
300 | since the DRS calibration is done by the C++ class CalFactFits
|
---|
301 |
|
---|
302 | perform the drs amplitude calibration of the event data
|
---|
303 | """
|
---|
304 | # shortcuts
|
---|
305 | blm = self.blm
|
---|
306 | gm = self.gm
|
---|
307 | tom = self.tom
|
---|
308 |
|
---|
309 | to_mV = 2000./4096.
|
---|
310 | #: 2D array with amplitude calibrated dat in mV
|
---|
311 | acal_data = self.data * to_mV # convert ADC counts to mV
|
---|
312 |
|
---|
313 |
|
---|
314 | for pixel in range( self.npix ):
|
---|
315 | #shortcuts
|
---|
316 | sc = self.start_cells[pixel]
|
---|
317 | roi = self.nroi
|
---|
318 | # rotate the pixel baseline mean to the Data startCell
|
---|
319 | acal_data[pixel,:] -= blm[pixel,sc:sc+roi]
|
---|
320 | # the 'trigger offset mean' does not need to be rolled
|
---|
321 | # on the contrary, it seems there is an offset in the DRS data,
|
---|
322 | # which is related to its distance to the startCell, not to its
|
---|
323 | # distance to the beginning of the physical pipeline in the DRS chip
|
---|
324 | acal_data[pixel,:] -= tom[pixel,0:roi]
|
---|
325 | # rotate the pixel gain mean to the Data startCell
|
---|
326 | acal_data[pixel,:] /= gm[pixel,sc:sc+roi]
|
---|
327 |
|
---|
328 |
|
---|
329 | self.acal_data = acal_data * 1907.35
|
---|
330 |
|
---|
331 | self.user_action_calib( self.acal_data,
|
---|
332 | np.reshape(self.data, (self.npix, self.nroi) ), blm, tom, gm, self.start_cells, self.nroi)
|
---|
333 |
|
---|
334 |
|
---|
335 | def baseline_read_values(self, file, bsl_hist='bsl_sum/hplt_mean'):
|
---|
336 | """
|
---|
337 | open ROOT file with baseline histogram and read baseline values
|
---|
338 |
|
---|
339 | *file* : name of the root file
|
---|
340 |
|
---|
341 | *bsl_hist* : path to the histogram containing the basline values
|
---|
342 | """
|
---|
343 |
|
---|
344 | try:
|
---|
345 | f = TFile(file)
|
---|
346 | except:
|
---|
347 | print 'Baseline data file could not be read: ', file
|
---|
348 | return
|
---|
349 |
|
---|
350 | h = f.Get(bsl_hist)
|
---|
351 |
|
---|
352 | for i in range(self.npix):
|
---|
353 | self.v_bsl[i] = h.GetBinContent(i+1)
|
---|
354 |
|
---|
355 | f.Close()
|
---|
356 |
|
---|
357 | def baseline_correct(self):
|
---|
358 | """ subtract baseline from the data
|
---|
359 |
|
---|
360 | DN 08.06.2011: I didn't use this function at all so far... don't know how well it works.
|
---|
361 | """
|
---|
362 |
|
---|
363 | for pixel in range(self.npix):
|
---|
364 | self.acal_data[pixel,:] -= self.v_bsl[pixel]
|
---|
365 |
|
---|
366 | def info(self):
|
---|
367 | """ print run information
|
---|
368 |
|
---|
369 | not very well implemented ... we need more info here.
|
---|
370 | """
|
---|
371 | print 'data file: ', self.data_file_name
|
---|
372 | print 'calib file: ', self.calib_file_name
|
---|
373 | print '... we need more information printed here ... '
|
---|
374 |
|
---|
375 | # -----------------------------------------------------------------------------
|
---|
376 | class RawDataFake( object ):
|
---|
377 | """ raw data FAKE access similar to real RawData access
|
---|
378 |
|
---|
379 | DO NOT USE ... its not working
|
---|
380 | """
|
---|
381 |
|
---|
382 |
|
---|
383 | def __init__(self, data_file_name, calib_file_name,
|
---|
384 | user_action_calib=lambda acal_data, data, blm, tom, gm, scells, nroi: None,
|
---|
385 | baseline_file_name=''):
|
---|
386 | self.__module__='pyfact'
|
---|
387 |
|
---|
388 | self.nroi = 300
|
---|
389 | self.npix = 9
|
---|
390 | self.nevents = 1000
|
---|
391 |
|
---|
392 | self.simulator = None
|
---|
393 |
|
---|
394 | self.time = np.ones(1024) * 0.5
|
---|
395 |
|
---|
396 |
|
---|
397 | self.event_id = c_ulong(0)
|
---|
398 | self.trigger_type = c_ushort(4)
|
---|
399 | self.data = np.zeros( self.npix * self.nroi, np.int16 ).reshape(self.npix ,self.nroi)
|
---|
400 | self.start_cells = np.zeros( self.npix, np.int16 )
|
---|
401 | self.board_times = np.zeros( 40, np.int32 )
|
---|
402 | def __iter__(self):
|
---|
403 | """ iterator """
|
---|
404 | return self
|
---|
405 |
|
---|
406 | def next(self):
|
---|
407 | """ used by __iter__ """
|
---|
408 | self.event_id = c_ulong(self.event_id.value + 1)
|
---|
409 | self.board_times = self.board_times + 42
|
---|
410 |
|
---|
411 | if self.event_id.value >= self.nevents:
|
---|
412 | raise StopIteration
|
---|
413 | else:
|
---|
414 | self._make_event_data()
|
---|
415 |
|
---|
416 | return self.__dict__
|
---|
417 |
|
---|
418 | def _make_event_data(self):
|
---|
419 | sample_times = self.time.cumsum() - time[0]
|
---|
420 |
|
---|
421 | # random start cell
|
---|
422 | self.start_cells = np.ones( self.npix, np.int16 ) * np.random.randint(0,1024)
|
---|
423 |
|
---|
424 | starttime = self.start_cells[0]
|
---|
425 |
|
---|
426 | signal = self._std_sinus_simu(sample_times, starttime)
|
---|
427 |
|
---|
428 | data = np.vstack( (signal,signal) )
|
---|
429 | for i in range(8):
|
---|
430 | data = np.vstack( (data,signal) )
|
---|
431 |
|
---|
432 | self.data = data
|
---|
433 |
|
---|
434 | def _std_sinus_simu(self, times, starttime):
|
---|
435 | period = 10 # in ns
|
---|
436 |
|
---|
437 | # give a jitter on starttime
|
---|
438 | starttime = np.random.normal(startime, 0.05)
|
---|
439 |
|
---|
440 | phase = 0.0
|
---|
441 | signal = 10 * np.sin(times * 2*np.pi/period + starttime + phase)
|
---|
442 |
|
---|
443 | # add some noise
|
---|
444 | noise = np.random.normal(0.0, 0.5, signal.shape)
|
---|
445 | signal += noise
|
---|
446 | return signal
|
---|
447 |
|
---|
448 | def info(self):
|
---|
449 | """ print run information
|
---|
450 |
|
---|
451 | """
|
---|
452 |
|
---|
453 | print 'data file: ', data_file_name
|
---|
454 | print 'calib file: ', calib_file_name
|
---|
455 | print 'calibration file'
|
---|
456 | print 'N baseline_mean: ', self.Nblm
|
---|
457 | print 'N gain mean: ', self.Ngm
|
---|
458 | print 'N TriggeroffsetMean: ', self.Ntom
|
---|
459 |
|
---|
460 | # -----------------------------------------------------------------------------
|
---|
461 | import ctypes
|
---|
462 |
|
---|
463 | class SlowData( FactFits ):
|
---|
464 | """ -Fact SlowData File-
|
---|
465 |
|
---|
466 | A Python wrapper for the fits-class implemented in factfits.h
|
---|
467 | provides easy access to the fits file meta data.
|
---|
468 |
|
---|
469 | * dictionary of file metadata - self.meta
|
---|
470 | * dict of table metadata - self.columns
|
---|
471 | * variable table column access, thus possibly increased speed while looping
|
---|
472 | """
|
---|
473 | def __init__(self, path):
|
---|
474 | """ creates meta and columns dictionaries
|
---|
475 | """
|
---|
476 | self.path = path
|
---|
477 | self.__module__ = 'pyfact'
|
---|
478 | try:
|
---|
479 | FactFits.__init__(self,path)
|
---|
480 | except IOError:
|
---|
481 | print 'problem accessing data file: ', data_file_name
|
---|
482 | raise # stop ! no data
|
---|
483 |
|
---|
484 | self.meta = self._make_meta_dict()
|
---|
485 | self.columns = self._make_columns_dict()
|
---|
486 |
|
---|
487 | self._treat_meta_dict()
|
---|
488 |
|
---|
489 |
|
---|
490 | # list of columns, which are already registered
|
---|
491 | # see method register()
|
---|
492 | self._registered_cols = []
|
---|
493 | # dict of column data, this is used, in order to be able to remove
|
---|
494 | # the ctypes of
|
---|
495 | self._table_cols = {}
|
---|
496 |
|
---|
497 | # I need to count the rows, since the normal loop mechanism seems not to work.
|
---|
498 | self._current_row = 0
|
---|
499 |
|
---|
500 | self.stacked_cols = {}
|
---|
501 |
|
---|
502 | def _make_meta_dict(self):
|
---|
503 | """ This method retrieves meta information about the fits file and
|
---|
504 | stores this information in a dict
|
---|
505 | return: dict
|
---|
506 | key: string - all capital letters
|
---|
507 | value: tuple( numerical value, string comment)
|
---|
508 | """
|
---|
509 | # intermediate variables for file metadata dict generation
|
---|
510 | keys=self.GetPy_KeyKeys()
|
---|
511 | values=self.GetPy_KeyValues()
|
---|
512 | comments=self.GetPy_KeyComments()
|
---|
513 | types=self.GetPy_KeyTypes()
|
---|
514 |
|
---|
515 | if len(keys) != len(values):
|
---|
516 | raise TypeError('len(keys)',len(keys),' != len(values)', len(values))
|
---|
517 | if len(keys) != len(types):
|
---|
518 | raise TypeError('len(keys)',len(keys),' != len(types)', len(types))
|
---|
519 | if len(keys) != len(comments):
|
---|
520 | raise TypeError('len(keys)',len(keys),' != len(comments)', len(comments))
|
---|
521 |
|
---|
522 | meta_dict = {}
|
---|
523 | for i in range(len(keys)):
|
---|
524 | type = types[i]
|
---|
525 | if type == 'I':
|
---|
526 | value = int(values[i])
|
---|
527 | elif type == 'F':
|
---|
528 | value = float(values[i])
|
---|
529 | elif type == 'B':
|
---|
530 | if values[i] == 'T':
|
---|
531 | value = True
|
---|
532 | elif values[i] == 'F':
|
---|
533 | value = False
|
---|
534 | else:
|
---|
535 | raise TypeError("meta-type is 'B', but meta-value is neither 'T' nor 'F'. meta-value:",values[i])
|
---|
536 | elif type == 'T':
|
---|
537 | value = values[i]
|
---|
538 | else:
|
---|
539 | raise TypeError("unknown meta-type: known meta types are: I,F,B and T. meta-type:",type)
|
---|
540 | meta_dict[keys[i]]=(value, comments[i])
|
---|
541 | return meta_dict
|
---|
542 |
|
---|
543 |
|
---|
544 | def _make_columns_dict(self):
|
---|
545 | """ This method retrieves information about the columns
|
---|
546 | stored inside the fits files internal binary table.
|
---|
547 | returns: dict
|
---|
548 | key: string column name -- all capital letters
|
---|
549 | values: tuple(
|
---|
550 | number of elements in table field - integer
|
---|
551 | size of element in bytes -- this is not really interesting for any user
|
---|
552 | might be ommited in future versions
|
---|
553 | type - a single character code -- should be translated into
|
---|
554 | a comrehensible word
|
---|
555 | unit - string like 'mV' or 'ADC count'
|
---|
556 | """
|
---|
557 | # intermediate variables for file table-metadata dict generation
|
---|
558 | keys=self.GetPy_ColumnKeys()
|
---|
559 | #offsets=self.GetPy_ColumnOffsets() #not needed on python level...
|
---|
560 | nums=self.GetPy_ColumnNums()
|
---|
561 | sizes=self.GetPy_ColumnSizes()
|
---|
562 | types=self.GetPy_ColumnTypes()
|
---|
563 | units=self.GetPy_ColumnUnits()
|
---|
564 |
|
---|
565 | # zip the values
|
---|
566 | values = zip(nums,sizes,types,units)
|
---|
567 | # create the columns dictionary
|
---|
568 | columns = dict(zip(keys ,values))
|
---|
569 | return columns
|
---|
570 |
|
---|
571 | def stack(self, on=True):
|
---|
572 | self.next()
|
---|
573 | for col in self._registered_cols:
|
---|
574 | if isinstance( self.dict[col], type(np.array('')) ):
|
---|
575 | self.stacked_cols[col] = self.dict[col]
|
---|
576 | else:
|
---|
577 | # elif isinstance(self.dict[col], ctypes._SimpleCData):
|
---|
578 | self.stacked_cols[col] = np.array(self.dict[col])
|
---|
579 | # else:
|
---|
580 | # raise TypeError("I don't know how to stack "+col+". It is of type: "+str(type(self.dict[col])))
|
---|
581 |
|
---|
582 | def register(self, col_name):
|
---|
583 | """ register for a column in the fits file
|
---|
584 |
|
---|
585 | after the call, this SlowData object will have a new member variable
|
---|
586 | self.col_name, if col_name is a key in self.colums
|
---|
587 |
|
---|
588 | the value will be updated after each call of next(), or while iterating over self.
|
---|
589 | NB: the initial value is zero(s)
|
---|
590 |
|
---|
591 | *col_name* : name of a key in self.columns, or 'all' to choose all.
|
---|
592 | """
|
---|
593 | columns = self.columns
|
---|
594 | if col_name.lower() == 'all':
|
---|
595 | for col in columns:
|
---|
596 | self._register(col)
|
---|
597 | else:
|
---|
598 | #check if colname is in columns:
|
---|
599 | if col_name not in columns:
|
---|
600 | error_msg = 'colname:'+ col_name +' is not a column in the binary table.\n'
|
---|
601 | error_msg+= 'possible colnames are\n'
|
---|
602 | for key in columns:
|
---|
603 | error_msg += key+' '
|
---|
604 | raise KeyError(error_msg)
|
---|
605 | else:
|
---|
606 | self._register(col_name)
|
---|
607 |
|
---|
608 | # 'private' method, do not use
|
---|
609 | def _register( self, colname):
|
---|
610 | columns = self.columns
|
---|
611 | local = None
|
---|
612 |
|
---|
613 | number_of_elements = int(columns[colname][0])
|
---|
614 | size_of_elements_in_bytes = int(columns[colname][1])
|
---|
615 | ctypecode_of_elements = columns[colname][2]
|
---|
616 | physical_unit_of_elements = columns[colname][3]
|
---|
617 |
|
---|
618 | # snippet from the C++ source code, or header file to be precise:
|
---|
619 | #case 'L': gLog << "bool(8)"; break;
|
---|
620 | #case 'B': gLog << "byte(8)"; break;
|
---|
621 | #case 'I': gLog << "short(16)"; break;
|
---|
622 | #case 'J': gLog << "int(32)"; break;
|
---|
623 | #case 'K': gLog << "int(64)"; break;
|
---|
624 | #case 'E': gLog << "float(32)"; break;
|
---|
625 | #case 'D': gLog << "double(64)"; break;
|
---|
626 |
|
---|
627 |
|
---|
628 |
|
---|
629 | # the fields inside the columns can either contain single numbers,
|
---|
630 | # or whole arrays of numbers as well.
|
---|
631 | # we treat single elements differently...
|
---|
632 | if number_of_elements == 1:
|
---|
633 | # allocate some memory for a single number according to its type
|
---|
634 | if ctypecode_of_elements == 'J': # J is for a 4byte int, i.e. an unsigned long
|
---|
635 | local = ctypes.c_ulong()
|
---|
636 | un_c_type = long
|
---|
637 | elif ctypecode_of_elements == 'I': # I is for a 2byte int, i.e. an unsinged int
|
---|
638 | local = ctypes.c_ushort()
|
---|
639 | un_c_type = int
|
---|
640 | elif ctypecode_of_elements == 'B': # B is for a byte
|
---|
641 | local = ctypes.c_ubyte()
|
---|
642 | un_c_type = int
|
---|
643 | elif ctypecode_of_elements == 'D':
|
---|
644 | local = ctypes.c_double()
|
---|
645 | un_c_type = float
|
---|
646 | elif ctypecode_of_elements == 'E':
|
---|
647 | local = ctypes.c_float()
|
---|
648 | un_c_type = float
|
---|
649 | elif ctypecode_of_elements == 'A':
|
---|
650 | local = ctypes.c_uchar()
|
---|
651 | un_c_type = chr
|
---|
652 | elif ctypecode_of_elements == 'K':
|
---|
653 | local = ctypes.c_ulonglong()
|
---|
654 | un_c_type = long
|
---|
655 | else:
|
---|
656 | raise TypeError('unknown ctypecode_of_elements:',ctypecode_of_elements)
|
---|
657 | else:
|
---|
658 | if ctypecode_of_elements == 'B': # B is for a byte
|
---|
659 | nptype = np.int8
|
---|
660 | elif ctypecode_of_elements == 'A': # A is for a char .. but I don't know how to handle it
|
---|
661 | nptype = np.int8
|
---|
662 | elif ctypecode_of_elements == 'I': # I is for a 2byte int
|
---|
663 | nptype = np.int16
|
---|
664 | elif ctypecode_of_elements == 'J': # J is for a 4byte int
|
---|
665 | nptype = np.int32
|
---|
666 | elif ctypecode_of_elements == 'K': # B is for a byte
|
---|
667 | nptype = np.int64
|
---|
668 | elif ctypecode_of_elements == 'E': # B is for a byte
|
---|
669 | nptype = np.float32
|
---|
670 | elif ctypecode_of_elements == 'D': # B is for a byte
|
---|
671 | nptype = np.float64
|
---|
672 | else:
|
---|
673 | raise TypeError('unknown ctypecode_of_elements:',ctypecode_of_elements)
|
---|
674 | local = np.zeros( number_of_elements, nptype)
|
---|
675 |
|
---|
676 | # Set the Pointer Address
|
---|
677 | self.SetPtrAddress(colname, local)
|
---|
678 | self._table_cols[colname] = local
|
---|
679 | if number_of_elements > 1:
|
---|
680 | self.__dict__[colname] = local
|
---|
681 | self.dict[colname] = local
|
---|
682 | else:
|
---|
683 | # remove any traces of ctypes:
|
---|
684 | self.__dict__[colname] = local.value
|
---|
685 | self.dict[colname] = local.value
|
---|
686 | self._registered_cols.append(colname)
|
---|
687 |
|
---|
688 |
|
---|
689 | def _treat_meta_dict(self):
|
---|
690 | """make 'interesting' meta information available like normal members.
|
---|
691 | non interesting are:
|
---|
692 | TFORM, TUNIT, and TTYPE
|
---|
693 | since these are available via the columns dict.
|
---|
694 | """
|
---|
695 |
|
---|
696 | self.number_of_rows = self.meta['NAXIS2'][0]
|
---|
697 | self.number_of_columns = self.meta['TFIELDS'][0]
|
---|
698 |
|
---|
699 | # there are some information in the meta dict, which are alsways there:
|
---|
700 | # there are regarded as not interesting:
|
---|
701 | uninteresting_meta = {}
|
---|
702 | uninteresting_meta['arraylike'] = {}
|
---|
703 | uninteresting = ['NAXIS', 'NAXIS1', 'NAXIS2',
|
---|
704 | 'TFIELDS',
|
---|
705 | 'XTENSION','EXTNAME','EXTREL',
|
---|
706 | 'BITPIX', 'PCOUNT', 'GCOUNT',
|
---|
707 | 'ORIGIN',
|
---|
708 | 'PACKAGE', 'COMPILED', 'CREATOR',
|
---|
709 | 'TELESCOP','TIMESYS','TIMEUNIT','VERSION']
|
---|
710 | for key in uninteresting:
|
---|
711 | if key in self.meta:
|
---|
712 | uninteresting_meta[key]=self.meta[key]
|
---|
713 | del self.meta[key]
|
---|
714 |
|
---|
715 | # the table meta data contains
|
---|
716 |
|
---|
717 |
|
---|
718 | # shortcut to access the meta dict. But this needs to
|
---|
719 | # be cleaned up quickly!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
---|
720 | meta = self.meta
|
---|
721 |
|
---|
722 | # loop over keys:
|
---|
723 | # * try to find array-like keys
|
---|
724 | arraylike = {}
|
---|
725 | singlelike = []
|
---|
726 | for key in self.meta:
|
---|
727 | stripped = key.rstrip('1234567890')
|
---|
728 | if stripped == key:
|
---|
729 | singlelike.append(key)
|
---|
730 | else:
|
---|
731 | if stripped not in arraylike:
|
---|
732 | arraylike[stripped] = 0
|
---|
733 | else:
|
---|
734 | arraylike[stripped] += 1
|
---|
735 | newmeta = {}
|
---|
736 | for key in singlelike:
|
---|
737 | newmeta[key.lower()] = meta[key]
|
---|
738 | for key in arraylike:
|
---|
739 | uninteresting_meta['arraylike'][key.lower()] = []
|
---|
740 | for i in range(arraylike[key]+1):
|
---|
741 | if key+str(i) in meta:
|
---|
742 | uninteresting_meta['arraylike'][key.lower()].append(meta[key+str(i)])
|
---|
743 | self.ui_meta = uninteresting_meta
|
---|
744 | # make newmeta self
|
---|
745 | for key in newmeta:
|
---|
746 | self.__dict__[key]=newmeta[key]
|
---|
747 |
|
---|
748 | dict = self.__dict__.copy()
|
---|
749 | del dict['meta']
|
---|
750 | del dict['ui_meta']
|
---|
751 | self.dict = dict
|
---|
752 |
|
---|
753 | def __iter__(self):
|
---|
754 | """ iterator """
|
---|
755 | return self
|
---|
756 |
|
---|
757 | def next(self):
|
---|
758 | """ use to iterate over the file
|
---|
759 |
|
---|
760 | do not forget to call register() before iterating over the file
|
---|
761 | call show() in order to find out, what parameters register() accepts.
|
---|
762 | or just call register('all') in case you are unsure.
|
---|
763 |
|
---|
764 | returns self
|
---|
765 | """
|
---|
766 | # Here one might check, if looping makes any sense, and if not
|
---|
767 | # one could stop looping or so...
|
---|
768 | # like this:
|
---|
769 | #
|
---|
770 | # if len(self._registered_cols) == 0:
|
---|
771 | # print 'warning: looping without any registered columns'
|
---|
772 | if self._current_row < self.number_of_rows:
|
---|
773 | if self.GetNextRow() == False:
|
---|
774 | raise StopIteration
|
---|
775 | for col in self._registered_cols:
|
---|
776 | if isinstance(self._table_cols[col], ctypes._SimpleCData):
|
---|
777 | self.__dict__[col] = self._table_cols[col].value
|
---|
778 | self.dict[col] = self._table_cols[col].value
|
---|
779 |
|
---|
780 | for col in self.stacked_cols:
|
---|
781 | if isinstance(self.dict[col], type(np.array(''))):
|
---|
782 | self.stacked_cols[col] = np.vstack( (self.stacked_cols[col],self.dict[col]) )
|
---|
783 | else:
|
---|
784 | self.stacked_cols[col] = np.vstack( (self.stacked_cols[col],np.array(self.dict[col])) )
|
---|
785 | self._current_row += 1
|
---|
786 | else:
|
---|
787 | raise StopIteration
|
---|
788 | return self
|
---|
789 |
|
---|
790 | def show(self):
|
---|
791 | """
|
---|
792 | """
|
---|
793 | pprint.pprint(self.dict)
|
---|
794 |
|
---|
795 |
|
---|
796 |
|
---|
797 |
|
---|
798 | class fnames( object ):
|
---|
799 | """ organize file names of a FACT data run
|
---|
800 |
|
---|
801 | """
|
---|
802 |
|
---|
803 | def __init__(self, specifier = ['012', '023', '2011', '11', '24'],
|
---|
804 | rpath = '/scratch_nfs/res/bsl/',
|
---|
805 | zipped = True):
|
---|
806 | """
|
---|
807 | specifier : list of strings defined as:
|
---|
808 | [ 'DRS calibration file', 'Data file', 'YYYY', 'MM', 'DD']
|
---|
809 |
|
---|
810 | rpath : directory path for the results; YYYYMMDD will be appended to rpath
|
---|
811 | zipped : use zipped (True) or unzipped (Data)
|
---|
812 |
|
---|
813 | """
|
---|
814 |
|
---|
815 | self.specifier = specifier
|
---|
816 | self.rpath = rpath
|
---|
817 | self.zipped = zipped
|
---|
818 |
|
---|
819 | self.make( self.specifier, self.rpath, self.zipped )
|
---|
820 |
|
---|
821 |
|
---|
822 | def make( self, specifier, rpath, zipped ):
|
---|
823 | """ create (make) the filenames
|
---|
824 |
|
---|
825 | names : dictionary of filenames, tags { 'data', 'drscal', 'results' }
|
---|
826 | data : name of the data file
|
---|
827 | drscal : name of the drs calibration file
|
---|
828 | results : radikal of file name(s) for results (to be completed by suffixes)
|
---|
829 | """
|
---|
830 |
|
---|
831 | self.specifier = specifier
|
---|
832 |
|
---|
833 | if zipped:
|
---|
834 | dpath = '/data00/fact-construction/raw/'
|
---|
835 | ext = '.fits.gz'
|
---|
836 | else:
|
---|
837 | dpath = '/data03/fact-construction/raw/'
|
---|
838 | ext = '.fits'
|
---|
839 |
|
---|
840 | year = specifier[2]
|
---|
841 | month = specifier[3]
|
---|
842 | day = specifier[4]
|
---|
843 |
|
---|
844 | yyyymmdd = year + month + day
|
---|
845 | dfile = specifier[1]
|
---|
846 | cfile = specifier[0]
|
---|
847 |
|
---|
848 | rpath = rpath + yyyymmdd + '/'
|
---|
849 | self.rpath = rpath
|
---|
850 | self.names = {}
|
---|
851 |
|
---|
852 | tmp = dpath + year + '/' + month + '/' + day + '/' + yyyymmdd + '_'
|
---|
853 | self.names['data'] = tmp + dfile + ext
|
---|
854 | self.names['drscal'] = tmp + cfile + '.drs' + ext
|
---|
855 | self.names['results'] = rpath + yyyymmdd + '_' + dfile + '_' + cfile
|
---|
856 |
|
---|
857 | self.data = self.names['data']
|
---|
858 | self.drscal = self.names['drscal']
|
---|
859 | self.results = self.names['results']
|
---|
860 |
|
---|
861 | def info( self ):
|
---|
862 | """ print complete filenames
|
---|
863 |
|
---|
864 | """
|
---|
865 |
|
---|
866 | print 'file names:'
|
---|
867 | print 'data: ', self.names['data']
|
---|
868 | print 'drs-cal: ', self.names['drscal']
|
---|
869 | print 'results: ', self.names['results']
|
---|
870 |
|
---|
871 | # end of class definition: fnames( object )
|
---|
872 |
|
---|
873 | def _test_SlowData( filename ):
|
---|
874 | print '-'*70
|
---|
875 | print "opened :", filename, " as 'file'"
|
---|
876 | print
|
---|
877 | print '-'*70
|
---|
878 | print 'type file.show() to look at its contents'
|
---|
879 | print "type file.register( columnname ) or file.register('all') in order to register columns"
|
---|
880 | print
|
---|
881 | print " due column-registration you declare, that you would like to retrieve the contents of one of the columns"
|
---|
882 | print " after column-registration, the 'file' has new member variables, they are named like the columns"
|
---|
883 | print " PLEASE NOTE: immediatly after registration, the members exist, but they are empty."
|
---|
884 | print " the values are assigned only, when you call file.next() or when you loop over the 'file'"
|
---|
885 | print
|
---|
886 | print "in order to loop over it, just go like this:"
|
---|
887 | print "for row in file:"
|
---|
888 | print " print row.columnname_one, row.columnname_two"
|
---|
889 | print
|
---|
890 | print ""
|
---|
891 | print '-'*70
|
---|
892 |
|
---|
893 |
|
---|
894 |
|
---|
895 | def _test_iter( nevents ):
|
---|
896 | """ test for function __iter__ """
|
---|
897 |
|
---|
898 | data_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_117.fits.gz'
|
---|
899 | calib_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_114.drs.fits.gz'
|
---|
900 | # data_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_028.fits.gz'
|
---|
901 | # calib_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_022.drs.fits.gz'
|
---|
902 | run = RawData( data_file_name, calib_file_name , return_dict=True)
|
---|
903 |
|
---|
904 | for event in run:
|
---|
905 | print 'ev ', event['event_id'], 'data[0,0] = ', event['acal_data'][0,0], 'start_cell[0] = ', event['start_cells'][0], 'trigger type = ', event['trigger_type']
|
---|
906 | if run.event_id == nevents:
|
---|
907 | break
|
---|
908 |
|
---|
909 | if __name__ == '__main__':
|
---|
910 | """ tests """
|
---|
911 | import sys
|
---|
912 | if len(sys.argv) == 1:
|
---|
913 | print 'showing test of iterator of RawData class'
|
---|
914 | print 'in order to test the SlowData classe please use:', sys.argv[0], 'fits-file-name'
|
---|
915 | _test_iter(10)
|
---|
916 |
|
---|
917 |
|
---|
918 | else:
|
---|
919 | print 'showing test of SlowData class'
|
---|
920 | print 'in case you wanted to test the RawData class, please give no commandline arguments'
|
---|
921 | file = SlowData(sys.argv[1])
|
---|
922 | _test_SlowData(sys.argv[1])
|
---|