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 | from scipy import signal
|
---|
9 |
|
---|
10 | # get the ROOT stuff + my shared libs
|
---|
11 | from 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
|
---|
15 | gSystem.Load('pyfits_h.so')
|
---|
16 | from ROOT import *
|
---|
17 |
|
---|
18 | class 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 |
|
---|
69 | class 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 | # -----------------------------------------------------------------------------
|
---|
299 | class 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 | # -----------------------------------------------------------------------------
|
---|
382 | class 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 |
|
---|
457 | def _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 |
|
---|
471 | if __name__ == '__main__':
|
---|
472 | """ tests """
|
---|
473 |
|
---|
474 | _test_iter(10)
|
---|