source: fact/tools/pyscripts/pyfact/pyfact.py@ 17730

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