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

Last change on this file since 13447 was 13441, checked in by neise, 13 years ago
init of class RawDataFake ... not yet functional
  • Property svn:executable set to *
File size: 16.5 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
8from scipy import signal
9
10# get the ROOT stuff + my shared libs
11from ROOT import gSystem
12# pyfits_h.so is made from pyfits.h and is used to access the data
13# make sure the location of pyfits_h.so is in LD_LIBRARY_PATH.
14# having it in PYTHONPATH is *not* sufficient
15gSystem.Load('pyfits_h.so')
16from ROOT import *
17
18class RawDataFeeder( object ):
19 """ Wrapper class for RawData class
20 capable of iterating over multiple RawData Files
21 """
22
23 def __init__(self, filelist):
24 """ *filelist* list of files to iterate over
25 the list should contain tuples, or sublists of two filenames
26 the first should be a data file (\*.fits.gz)
27 the second should be an amplitude calibration file(\*.drs.fits.gz)
28 """
29 # sanity check for input
30 if type(filelist) != type(list()):
31 raise TypeError('filelist should be a list')
32 for entry in filelist:
33 if len(entry) != 2:
34 raise TypeError('the entries of filelist should have length == 2')
35 for path in entry:
36 if type(path) != type(str()):
37 raise TypeError('the entries of filelist should be path, i.e. of type str()')
38 #todo check if 'path' is a valid path
39 # else: throw an Exception, or Warning?
40
41 self.filelist = filelist
42 self._current_RawData = RawData(filelist[0][0], filelist[0][1], return_dict=True)
43 del filelist[0]
44
45 def __iter__(self):
46 return self
47
48 def next():
49 """ Method being called by the iterator.
50 Since the RawData Objects are simply looped over, the event_id from the
51 RawData object will not be unique.
52 Each RawData obejct will start with event_id = 1 as usual.
53 """
54 try:
55 return self._current_RawData.next()
56 except StopIteration:
57 # current_RawData was completely processed
58 # delete it (I hope this calls the destructor of the fits file and/or closes it)
59 del self._current_RawData
60 # and remake it, if possible
61 if len(self.filelist) > 0:
62 self._current_RawData = RawData(filelist[0][0], filelist[0][1], return_dict=True)
63 del filelist[0]
64 else:
65 raise
66
67
68
69class RawData( object ):
70 """ raw data access and calibration
71
72 - open raw data file and drs calibration file
73 - performs amplitude calibration
74 - performs baseline substraction if wanted
75 - provides all data in an array:
76 row = number of pixel
77 col = length of region of interest
78
79 """
80
81
82 def __init__(self, data_file_name, calib_file_name,
83 user_action_calib=lambda acal_data, data, blm, tom, gm, scells, nroi: None,
84 baseline_file_name='',
85 return_dict = None,
86 do_calibration = True):
87 """ initialize object
88
89 open data file and calibration data file
90 get basic information about the data in data_file_name
91 allocate buffers for data access
92
93 data_file_name : fits or fits.gz file of the data including the path
94 calib_file_name : fits or fits.gz file containing DRS calibration data
95 baseline_file_name : npy file containing the baseline values
96 """
97 self.__module__='pyfact'
98 # manual implementation of default value, but I need to find out
99 # if the user of this class is aware of the new option
100 if return_dict == None:
101 print 'Warning: Rawdata.__init__() has a new option "return_dict"'
102 print 'the default value of this option is False, so nothing changes for you at the moment.'
103 print
104 print 'you probably want, to get a dictionary out of the next() method anyway'
105 print ' so please change your scripts and set this option to True, for the moment'
106 print 'e.g. like this: run = RawData(data_filename, calib_filename, return_dict = True)'
107 print "after a while, the default value, will turn to True .. so you don't have to give the option anymore"
108 print 'and some time later, the option will not be supported anymore'
109 return_dict = False
110 self.return_dict = return_dict
111
112 self.do_calibration = do_calibration
113
114 self.data_file_name = data_file_name
115 self.calib_file_name = calib_file_name
116 self.baseline_file_name = baseline_file_name
117
118 self.user_action_calib = user_action_calib
119
120 # baseline correction: True / False
121 if len(baseline_file_name) == 0:
122 self.correct_baseline = False
123 else:
124 self.correct_baseline = True
125
126 # access data file
127 try:
128 data_file = fits(self.data_file_name)
129 except IOError:
130 print 'problem accessing data file: ', data_file_name
131 raise # stop ! no data
132 #: data file (fits object)
133 self.data_file = data_file
134
135 # get basic information about the data file
136 #: region of interest (number of DRS slices read)
137 self.nroi = data_file.GetUInt('NROI')
138 #: number of pixels (should be 1440)
139 self.npix = data_file.GetUInt('NPIX')
140 #: number of events in the data run
141 self.nevents = data_file.GetNumRows()
142
143 # allocate the data memories
144 self.event_id = c_ulong()
145 self.trigger_type = c_ushort()
146 #: 1D array with raw data
147 self.data = np.zeros( self.npix * self.nroi, np.int16 ).reshape(self.npix ,self.nroi)
148 #: slice where drs readout started
149 self.start_cells = np.zeros( self.npix, np.int16 )
150 #: time when the FAD was triggered, in some strange units...
151 self.board_times = np.zeros( 40, np.int32 )
152
153 # set the pointers to the data++
154 data_file.SetPtrAddress('EventNum', self.event_id)
155 data_file.SetPtrAddress('TriggerType', self.trigger_type)
156 data_file.SetPtrAddress('StartCellData', self.start_cells)
157 data_file.SetPtrAddress('Data', self.data)
158 data_file.SetPtrAddress('BoardTime', self.board_times)
159
160 # open the calibration file
161 try:
162 calib_file = fits(self.calib_file_name)
163 except IOError:
164 print 'problem accessing calibration file: ', calib_file_name
165 raise
166 #: drs calibration file
167 self.calib_file = calib_file
168
169 baseline_mean = calib_file.GetN('BaselineMean')
170 gain_mean = calib_file.GetN('GainMean')
171 trigger_offset_mean = calib_file.GetN('TriggerOffsetMean')
172
173 self.Nblm = baseline_mean / self.npix
174 self.Ngm = gain_mean / self.npix
175 self.Ntom = trigger_offset_mean / self.npix
176
177 self.blm = np.zeros(baseline_mean, np.float32).reshape(self.npix , self.Nblm)
178 self.gm = np.zeros(gain_mean, np.float32).reshape(self.npix , self.Ngm)
179 self.tom = np.zeros(trigger_offset_mean, np.float32).reshape(self.npix , self.Ntom)
180
181 calib_file.SetPtrAddress('BaselineMean', self.blm)
182 calib_file.SetPtrAddress('GainMean', self.gm)
183 calib_file.SetPtrAddress('TriggerOffsetMean', self.tom)
184 calib_file.GetRow(0)
185
186 # make calibration constants double, so we never need to roll
187 self.blm = np.hstack((self.blm, self.blm))
188 self.gm = np.hstack((self.gm, self.gm))
189 self.tom = np.hstack((self.tom, self.tom))
190
191 self.v_bsl = np.zeros(self.npix) # array of baseline values (all ZERO)
192
193 def __iter__(self):
194 """ iterator """
195 return self
196
197 def __add__(self, jump_over):
198 self.data_file.GetRow(jump_over)
199 return self
200
201 def next(self):
202 """ used by __iter__ """
203 if self.data_file.GetNextRow() == False:
204 raise StopIteration
205 else:
206 if self.do_calibration == True:
207 self.calibrate_drs_amplitude()
208
209 #print 'nevents = ', self.nevents, 'event_id = ', self.event_id.value
210 if self.return_dict:
211 return self.__dict__
212 else:
213 return self.acal_data, self.start_cells, self.trigger_type.value
214
215 def next_event(self):
216 """ load the next event from disk and calibrate it
217 """
218 self.data_file.GetNextRow()
219 self.calibrate_drs_amplitude()
220
221 def calibrate_drs_amplitude(self):
222 """ perform the drs amplitude calibration of the event data
223
224 """
225 # shortcuts
226 blm = self.blm
227 gm = self.gm
228 tom = self.tom
229
230 to_mV = 2000./4096.
231 #: 2D array with amplitude calibrated dat in mV
232 acal_data = self.data * to_mV # convert ADC counts to mV
233
234
235 for pixel in range( self.npix ):
236 #shortcuts
237 sc = self.start_cells[pixel]
238 roi = self.nroi
239 # rotate the pixel baseline mean to the Data startCell
240 acal_data[pixel,:] -= blm[pixel,sc:sc+roi]
241 # the 'trigger offset mean' does not need to be rolled
242 # on the contrary, it seems there is an offset in the DRS data,
243 # which is related to its distance to the startCell, not to its
244 # distance to the beginning of the physical pipeline in the DRS chip
245 acal_data[pixel,:] -= tom[pixel,0:roi]
246 # rotate the pixel gain mean to the Data startCell
247 acal_data[pixel,:] /= gm[pixel,sc:sc+roi]
248
249
250 self.acal_data = acal_data * 1907.35
251
252 self.user_action_calib( self.acal_data,
253 np.reshape(self.data, (self.npix, self.nroi) ), blm, tom, gm, self.start_cells, self.nroi)
254
255
256 def baseline_read_values(self, file, bsl_hist='bsl_sum/hplt_mean'):
257 """
258
259 open ROOT file with baseline histogram and read baseline values
260 file name of the root file
261 bsl_hist path to the histogram containing the basline values
262
263 """
264
265 try:
266 f = TFile(file)
267 except:
268 print 'Baseline data file could not be read: ', file
269 return
270
271 h = f.Get(bsl_hist)
272
273 for i in range(self.npix):
274 self.v_bsl[i] = h.GetBinContent(i+1)
275
276 f.Close()
277
278 def baseline_correct(self):
279 """ subtract baseline from the data
280
281 """
282
283 for pixel in range(self.npix):
284 self.acal_data[pixel,:] -= self.v_bsl[pixel]
285
286 def info(self):
287 """ print run information
288
289 """
290
291 print 'data file: ', data_file_name
292 print 'calib file: ', calib_file_name
293 print 'calibration file'
294 print 'N baseline_mean: ', self.Nblm
295 print 'N gain mean: ', self.Ngm
296 print 'N TriggeroffsetMean: ', self.Ntom
297
298# -----------------------------------------------------------------------------
299class RawDataFake( object ):
300 """ raw data FAKE access similar to real RawData access
301 """
302
303
304 def __init__(self, data_file_name, calib_file_name,
305 user_action_calib=lambda acal_data, data, blm, tom, gm, scells, nroi: None,
306 baseline_file_name=''):
307 self.__module__='pyfact'
308
309 self.nroi = 300
310 self.npix = 9
311 self.nevents = 1000
312
313 self.simulator = None
314
315 self.time = np.ones(1024) * 0.5
316
317
318 self.event_id = c_ulong(0)
319 self.trigger_type = c_ushort(4)
320 self.data = np.zeros( self.npix * self.nroi, np.int16 ).reshape(self.npix ,self.nroi)
321 self.start_cells = np.zeros( self.npix, np.int16 )
322 self.board_times = np.zeros( 40, np.int32 )
323 def __iter__(self):
324 """ iterator """
325 return self
326
327 def next(self):
328 """ used by __iter__ """
329 self.event_id = c_ulong(self.event_id.value + 1)
330 self.board_times = self.board_times + 42
331
332 if self.event_id.value >= self.nevents:
333 raise StopIteration
334 else:
335 self._make_event_data()
336
337 return self.__dict__
338
339 def _make_event_data(self):
340 sample_times = self.time.cumsum() - time[0]
341
342 # random start cell
343 self.start_cells = np.ones( self.npix, np.int16 ) * np.random.randint(0,1024)
344
345 starttime = self.start_cells[0]
346
347 signal = self._std_sinus_simu(sample_times, starttime)
348
349 data = np.vstack( (signal,signal) )
350 for i in range(8):
351 data = np.vstack( (data,signal) )
352
353 self.data = data
354
355 def _std_sinus_simu(self, times, starttime):
356 period = 10 # in ns
357
358 # give a jitter on starttime
359 starttime = np.random.normal(startime, 0.05)
360
361 phase = 0.0
362 signal = 10 * np.sin(times * 2*np.pi/period + starttime + phase)
363
364 # add some noise
365 noise = np.random.normal(0.0, 0.5, signal.shape)
366 signal += noise
367 return signal
368
369 def info(self):
370 """ print run information
371
372 """
373
374 print 'data file: ', data_file_name
375 print 'calib file: ', calib_file_name
376 print 'calibration file'
377 print 'N baseline_mean: ', self.Nblm
378 print 'N gain mean: ', self.Ngm
379 print 'N TriggeroffsetMean: ', self.Ntom
380
381# -----------------------------------------------------------------------------
382class fnames( object ):
383 """ organize file names of a FACT data run
384
385 """
386
387 def __init__(self, specifier = ['012', '023', '2011', '11', '24'],
388 rpath = '/scratch_nfs/res/bsl/',
389 zipped = True):
390 """
391 specifier : list of strings defined as:
392 [ 'DRS calibration file', 'Data file', 'YYYY', 'MM', 'DD']
393
394 rpath : directory path for the results; YYYYMMDD will be appended to rpath
395 zipped : use zipped (True) or unzipped (Data)
396
397 """
398
399 self.specifier = specifier
400 self.rpath = rpath
401 self.zipped = zipped
402
403 self.make( self.specifier, self.rpath, self.zipped )
404
405
406 def make( self, specifier, rpath, zipped ):
407 """ create (make) the filenames
408
409 names : dictionary of filenames, tags { 'data', 'drscal', 'results' }
410 data : name of the data file
411 drscal : name of the drs calibration file
412 results : radikal of file name(s) for results (to be completed by suffixes)
413 """
414
415 self.specifier = specifier
416
417 if zipped:
418 dpath = '/data00/fact-construction/raw/'
419 ext = '.fits.gz'
420 else:
421 dpath = '/data03/fact-construction/raw/'
422 ext = '.fits'
423
424 year = specifier[2]
425 month = specifier[3]
426 day = specifier[4]
427
428 yyyymmdd = year + month + day
429 dfile = specifier[1]
430 cfile = specifier[0]
431
432 rpath = rpath + yyyymmdd + '/'
433 self.rpath = rpath
434 self.names = {}
435
436 tmp = dpath + year + '/' + month + '/' + day + '/' + yyyymmdd + '_'
437 self.names['data'] = tmp + dfile + ext
438 self.names['drscal'] = tmp + cfile + '.drs' + ext
439 self.names['results'] = rpath + yyyymmdd + '_' + dfile + '_' + cfile
440
441 self.data = self.names['data']
442 self.drscal = self.names['drscal']
443 self.results = self.names['results']
444
445 def info( self ):
446 """ print complete filenames
447
448 """
449
450 print 'file names:'
451 print 'data: ', self.names['data']
452 print 'drs-cal: ', self.names['drscal']
453 print 'results: ', self.names['results']
454
455# end of class definition: fnames( object )
456
457def _test_iter( nevents ):
458 """ test for function __iter__ """
459
460 data_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_117.fits.gz'
461 calib_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_114.drs.fits.gz'
462# data_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_028.fits.gz'
463# calib_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_022.drs.fits.gz'
464 run = RawData( data_file_name, calib_file_name , return_dict=True)
465
466 for event in run:
467 print 'ev ', event['event_id'].value, 'data[0,0] = ', event['acal_data'][0,0], 'start_cell[0] = ', event['start_cells'][0], 'trigger type = ', event['trigger_type']
468 if run.event_id.value == nevents:
469 break
470
471if __name__ == '__main__':
472 """ tests """
473
474 _test_iter(10)
Note: See TracBrowser for help on using the repository browser.