| 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 fnames( object ):
|
|---|
| 300 | """ organize file names of a FACT data run
|
|---|
| 301 |
|
|---|
| 302 | """
|
|---|
| 303 |
|
|---|
| 304 | def __init__(self, specifier = ['012', '023', '2011', '11', '24'],
|
|---|
| 305 | rpath = '/scratch_nfs/res/bsl/',
|
|---|
| 306 | zipped = True):
|
|---|
| 307 | """
|
|---|
| 308 | specifier : list of strings defined as:
|
|---|
| 309 | [ 'DRS calibration file', 'Data file', 'YYYY', 'MM', 'DD']
|
|---|
| 310 |
|
|---|
| 311 | rpath : directory path for the results; YYYYMMDD will be appended to rpath
|
|---|
| 312 | zipped : use zipped (True) or unzipped (Data)
|
|---|
| 313 |
|
|---|
| 314 | """
|
|---|
| 315 |
|
|---|
| 316 | self.specifier = specifier
|
|---|
| 317 | self.rpath = rpath
|
|---|
| 318 | self.zipped = zipped
|
|---|
| 319 |
|
|---|
| 320 | self.make( self.specifier, self.rpath, self.zipped )
|
|---|
| 321 |
|
|---|
| 322 |
|
|---|
| 323 | def make( self, specifier, rpath, zipped ):
|
|---|
| 324 | """ create (make) the filenames
|
|---|
| 325 |
|
|---|
| 326 | names : dictionary of filenames, tags { 'data', 'drscal', 'results' }
|
|---|
| 327 | data : name of the data file
|
|---|
| 328 | drscal : name of the drs calibration file
|
|---|
| 329 | results : radikal of file name(s) for results (to be completed by suffixes)
|
|---|
| 330 | """
|
|---|
| 331 |
|
|---|
| 332 | self.specifier = specifier
|
|---|
| 333 |
|
|---|
| 334 | if zipped:
|
|---|
| 335 | dpath = '/data00/fact-construction/raw/'
|
|---|
| 336 | ext = '.fits.gz'
|
|---|
| 337 | else:
|
|---|
| 338 | dpath = '/data03/fact-construction/raw/'
|
|---|
| 339 | ext = '.fits'
|
|---|
| 340 |
|
|---|
| 341 | year = specifier[2]
|
|---|
| 342 | month = specifier[3]
|
|---|
| 343 | day = specifier[4]
|
|---|
| 344 |
|
|---|
| 345 | yyyymmdd = year + month + day
|
|---|
| 346 | dfile = specifier[1]
|
|---|
| 347 | cfile = specifier[0]
|
|---|
| 348 |
|
|---|
| 349 | rpath = rpath + yyyymmdd + '/'
|
|---|
| 350 | self.rpath = rpath
|
|---|
| 351 | self.names = {}
|
|---|
| 352 |
|
|---|
| 353 | tmp = dpath + year + '/' + month + '/' + day + '/' + yyyymmdd + '_'
|
|---|
| 354 | self.names['data'] = tmp + dfile + ext
|
|---|
| 355 | self.names['drscal'] = tmp + cfile + '.drs' + ext
|
|---|
| 356 | self.names['results'] = rpath + yyyymmdd + '_' + dfile + '_' + cfile
|
|---|
| 357 |
|
|---|
| 358 | self.data = self.names['data']
|
|---|
| 359 | self.drscal = self.names['drscal']
|
|---|
| 360 | self.results = self.names['results']
|
|---|
| 361 |
|
|---|
| 362 | def info( self ):
|
|---|
| 363 | """ print complete filenames
|
|---|
| 364 |
|
|---|
| 365 | """
|
|---|
| 366 |
|
|---|
| 367 | print 'file names:'
|
|---|
| 368 | print 'data: ', self.names['data']
|
|---|
| 369 | print 'drs-cal: ', self.names['drscal']
|
|---|
| 370 | print 'results: ', self.names['results']
|
|---|
| 371 |
|
|---|
| 372 | # end of class definition: fnames( object )
|
|---|
| 373 |
|
|---|
| 374 | def _test_iter( nevents ):
|
|---|
| 375 | """ test for function __iter__ """
|
|---|
| 376 |
|
|---|
| 377 | data_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_117.fits.gz'
|
|---|
| 378 | calib_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_114.drs.fits.gz'
|
|---|
| 379 | # data_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_028.fits.gz'
|
|---|
| 380 | # calib_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_022.drs.fits.gz'
|
|---|
| 381 | run = RawData( data_file_name, calib_file_name , return_dict=True)
|
|---|
| 382 |
|
|---|
| 383 | for event in run:
|
|---|
| 384 | 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']
|
|---|
| 385 | if run.event_id.value == nevents:
|
|---|
| 386 | break
|
|---|
| 387 |
|
|---|
| 388 | if __name__ == '__main__':
|
|---|
| 389 | """ tests """
|
|---|
| 390 |
|
|---|
| 391 | _test_iter(10)
|
|---|
