Index: /fact/tools/pyscripts/pyfact/pyfact.py =================================================================== --- /fact/tools/pyscripts/pyfact/pyfact.py (revision 13123) +++ /fact/tools/pyscripts/pyfact/pyfact.py (revision 13124) @@ -4,5 +4,5 @@ # ETH Zurich # -from ctypes import * +from ctypes import * import numpy as np from scipy import signal @@ -14,7 +14,8 @@ from ROOT import * -class rawdata( object ): - """ - raw data access and calibration + +class RawData( object ): + """ raw data access and calibration + - open raw data file and drs calibration file - performs amplitude calibration @@ -23,22 +24,30 @@ row = number of pixel col = length of region of interest + """ - # constructor of the classe - def __init__( self, dfname, calfname, bslfname='' ): - """ + + def __init__(self, data_file_name, calib_file_name, + user_action_calib=lambda acal_data, data, blm, tom, gm, scells, nroi: None, + baseline_file_name=''): + """ initialize object + open data file and calibration data file - get basic information about the data in dfname + get basic information about the data in data_file_name allocate buffers for data access - dfname : fits or fits.gz file containing the data including the path - calfname : fits or fits.gz file containing DRS calibration data - bslfname : npy file containing the baseline values - """ - self.dfname = dfname - self.calfname = calfname - self.bslfname = bslfname + data_file_name : fits or fits.gz file of the data including the path + calib_file_name : fits or fits.gz file containing DRS calibration data + baseline_file_name : npy file containing the baseline values + + """ + + self.data_file_name = data_file_name + self.calib_file_name = calib_file_name + self.baseline_file_name = baseline_file_name + + self.user_action_calib = user_action_calib # baseline correction: True / False - if len( bslfname ) == 0: + if len(baseline_file_name) == 0: self.correct_baseline = False else: @@ -47,208 +56,175 @@ # access data file try: - df = fits( self.dfname ) + data_file = fits(self.data_file_name) except IOError: - print 'problem accessing data file: ', dfname - raise # stop ! no data - self.df = df + print 'problem accessing data file: ', data_file_name + raise # stop ! no data + #: data file (fits object) + self.data_file = data_file # get basic information about the data file - self.NROI = df.GetUInt( 'NROI' ) # region of interest (length of DRS pipeline read out) - self.NPIX = df.GetUInt( 'NPIX' ) # number of pixels (should be 1440) - self.NEvents = df.GetNumRows() # find number of events + #: region of interest (number of DRS slices read) + self.nroi = data_file.GetUInt('NROI') + #: number of pixels (should be 1440) + self.npix = data_file.GetUInt('NPIX') + #: number of events in the data run + self.nevents = data_file.GetNumRows() + # allocate the data memories - self.evNum = c_ulong() - self.trigType = c_ushort() - self.Data = np.zeros( self.NPIX * self.NROI, np.int16 ) - self.startCells = np.zeros( self.NPIX, np.int16 ) + self.event_id = c_ulong() + self.trigger_type = c_ushort() + #: 1D array with raw data + self.data = np.zeros( self.npix * self.nroi, np.int16 ) + #: slice where drs readout started + self.start_cells = np.zeros( self.npix, np.int16 ) + #: time when the FAD was triggered, in some strange units... + self.board_times = np.zeros( 40, np.int32 ) + # set the pointers to the data++ - df.SetPtrAddress( 'EventNum', self.evNum ) - df.SetPtrAddress( 'TriggerType', self.trigType ) - df.SetPtrAddress( 'StartCellData', self.startCells ) # DRS readout start cell - df.SetPtrAddress( 'Data', self.Data ) # this is what you would expect - # df.GetNextRow() # access the first event - - # access calibration file + data_file.SetPtrAddress('EventNum', self.event_id) + data_file.SetPtrAddress('TriggerType', self.trigger_type) + data_file.SetPtrAddress('StartCellData', self.start_cells) + data_file.SetPtrAddress('Data', self.data) + data_file.SetPtrAddress('BoardTime', self.board_times) + + # open the calibration file try: - calf = fits( self.calfname ) + calib_file = fits(self.calib_file_name) except IOError: - print 'problem accessing calibration file: ', calfname + print 'problem accessing calibration file: ', calib_file_name raise - self.calf = calf - # - BaselineMean = calf.GetN('BaselineMean') - GainMean = calf.GetN('GainMean') - TriggerOffsetMean = calf.GetN('TriggerOffsetMean') - - self.blm = np.zeros( BaselineMean, np.float32 ) - self.gm = np.zeros( GainMean, np.float32 ) - self.tom = np.zeros( TriggerOffsetMean, np.float32 ) - - self.Nblm = BaselineMean / self.NPIX - self.Ngm = GainMean / self.NPIX - self.Ntom = TriggerOffsetMean / self.NPIX - - calf.SetPtrAddress( 'BaselineMean', self.blm ) - calf.SetPtrAddress( 'GainMean', self.gm ) - calf.SetPtrAddress( 'TriggerOffsetMean', self.tom ) - calf.GetRow(0) - - self.v_bsl = np.zeros( self.NPIX ) # array with baseline values (all ZERO) - self.smoothData = None - self.maxPos = None - self.maxAmp = None - - - def next( self ): - """ - load the next event from disk and calibrate it - """ - self.df.GetNextRow() - self.calibrate_drsAmplitude() - - - def calibrate_drsAmplitude( self ): - """ - perform amplitude calibration for the event - """ - tomV = 2000./4096. - acalData = self.Data * tomV # convert into mV - - # reshape arrays: row = pixel, col = drs_slice - acalData = np.reshape( acalData, (self.NPIX, self.NROI) ) - blm = np.reshape( self.blm, (self.NPIX, 1024) ) - tom = np.reshape( self.tom, (self.NPIX, 1024) ) - gm = np.reshape( self.gm, (self.NPIX, 1024) ) - - # print 'acal Data ', acalData.shape - # print 'blm shape ', blm.shape - # print 'gm shape ', gm.shape - - for pixel in range( self.NPIX ): + #: drs calibration file + self.calib_file = calib_file + + baseline_mean = calib_file.GetN('BaselineMean') + gain_mean = calib_file.GetN('GainMean') + trigger_offset_mean = calib_file.GetN('TriggerOffsetMean') + + self.blm = np.zeros(baseline_mean, np.float32) + self.gm = np.zeros(gain_mean, np.float32) + self.tom = np.zeros(trigger_offset_mean, np.float32) + + self.Nblm = baseline_mean / self.npix + self.Ngm = gain_mean / self.npix + self.Ntom = trigger_offset_mean / self.npix + + calib_file.SetPtrAddress('BaselineMean', self.blm) + calib_file.SetPtrAddress('GainMean', self.gm) + calib_file.SetPtrAddress('TriggerOffsetMean', self.tom) + calib_file.GetRow(0) + + self.v_bsl = np.zeros(self.npix) # array of baseline values (all ZERO) + self.data_saverage_out = None + self.pulse_time_of_maximum = None + self.pulse_amplitude = None + + def __iter__(self): + """ iterator """ + return self + + def __add__(self, jump_over): + self.data_file.GetRow(jump_over) + return self + + def next(self): + """ used by __iter__ """ + + if self.data_file.GetNextRow() == False: + raise StopIteration + else: + self.calibrate_drs_amplitude() + + #print 'nevents = ', self.nevents, 'event_id = ', self.event_id.value + + return self.acal_data, self.start_cells, self.trigger_type.value + + def next_event(self): + """ load the next event from disk and calibrate it + + """ + + self.data_file.GetNextRow() + self.calibrate_drs_amplitude() + + def calibrate_drs_amplitude(self): + """ perform the drs amplitude calibration of the event data + + """ + + to_mV = 2000./4096. + #: 2D array with amplitude calibrated dat in mV + acal_data = self.data * to_mV # convert ADC counts to mV + + # make 2D arrays: row = pixel, col = drs_slice + acal_data = np.reshape(acal_data, (self.npix, self.nroi) ) + blm = np.reshape(self.blm, (self.npix, self.Nblm) ) + tom = np.reshape(self.tom, (self.npix, self.Ntom) ) + gm = np.reshape(self.gm, (self.npix, self.Ngm) ) + + for pixel in range( self.npix ): # rotate the pixel baseline mean to the Data startCell - blm_pixel = np.roll( blm[pixel,:], -self.startCells[pixel] ) - acalData[pixel,:] -= blm_pixel[0:self.NROI] - acalData[pixel,:] -= tom[pixel, 0:self.NROI] - acalData[pixel,:] /= gm[pixel, 0:self.NROI] + blm_pixel = np.roll( blm[pixel,:], -self.start_cells[pixel] ) + tom_pixel = np.roll( tom[pixel,:], -self.start_cells[pixel] ) + gm_pixel = np.roll( gm[pixel,:], -self.start_cells[pixel] ) + acal_data[pixel,:] -= blm_pixel[0:self.nroi] + acal_data[pixel,:] -= tom_pixel[0:self.nroi] + acal_data[pixel,:] /= gm_pixel[0:self.nroi] - self.acalData = acalData * 1907.35 - - # print 'acalData ', self.acalData[0:2,0:20] - - def filterSlidingAverage( self , windowSize = 4): - """ sliding average filter - using: - self.acalData - filling array: - self.smoothData - """ - #scipy.signal.lfilter(b, a, x, axis=-1, zi=None) - smoothData = self.acalData.copy() - b = np.ones( windowSize ) - a = np.zeros( windowSize ) - a[0] = len(b) - smoothData[:,:] = signal.lfilter(b, a, smoothData[:,:]) - - self.smoothData = smoothData - - def filterCFD( self, length=10, ratio=0.75): - """ constant fraction filter - using: - self.smoothData - filling array: - self.cfdData - """ - if self.smoothData == None: - print 'error pyfact.filterCFD was called without prior call to filterSlidingAverage' - print ' variable self.smoothData is needed ' - pass - cfdData = self.smoothData.copy() - b = np.zeros( length ) - a = np.zeros( length ) - b[0] = -1. * ratio - b[length-1] = 1. - a[0] = 1. - cfdData[:,:] = signal.lfilter(b, a, cfdData[:,:]) - - self.cfdData = cfdData - - def findPeak (self, min=30, max=250): - """ find maximum in search window - using: - self.smoothData - filling arrays: - self.maxPos - self.maxAmp - """ - if self.smoothData == None: - print 'error pyfact.findPeakMax was called without prior call to filterSlidingAverage' - print ' variable self.smoothData is needed ' - pass - maxPos = np.argmax( self.smoothData[:,min:max] , 1) - maxAmp = np.max( self.smoothData[:,min:max] , 1) - self.maxPos = maxPos - self.maxAmp = maxAmp - - def sumAroundPeak (self, left=13, right=23): - """ integrate signal in gate around Peak - using: - self.maxPos - self.acalData - filling array: - self.integral - """ - if self.maxPos == None: - print 'error pyfact.sumAroundPeak was called without prior call of findPeak' - print ' variable self.maxPos is needed' - pass - - sums = np.empty( self.NPIX ) - for pixel in range( self.NPIX ): - min = self.maxPos[pixel]-left - max = self.maxPos[pixel]+right - sums[pixel] = self.acalData[pixel,min:max].sum() - - self.integral = sums - - def ReadBaseline( self, file, bsl_hist = 'bsl_sum/hplt_mean' ): - """ + self.acal_data = acal_data * 1907.35 + + #print 'blm _pyfact', blm[0,0:20] + #t = np.roll( blm[0,:], -self.start_cells[0] ) + #print 'blm _pyfact', t[0:20] + #print 'start_pyfact: ', self.start_cells[0] + #print 'acal _pyfact: ', self.acal_data[0,0:10] + #t = np.roll( gm[0,:], -self.start_cells[0] ) + #print 'gm _pyfact: ', t[0:10] + self.user_action_calib( self.acal_data, + np.reshape(self.data, (self.npix, self.nroi) ), blm, tom, gm, self.start_cells, self.nroi) + + + def baseline_read_values(self, file, bsl_hist='bsl_sum/hplt_mean'): + """ + open ROOT file with baseline histogram and read baseline values file name of the root file bsl_hist path to the histogram containing the basline values - """ + + """ + try: - f = TFile( file ) + f = TFile(file) except: print 'Baseline data file could not be read: ', file return - h = f.Get( bsl_hist ) - - for i in range( self.NPIX ): - self.v_bsl[i] = h.GetBinContent( i+1 ) + h = f.Get(bsl_hist) + + for i in range(self.npix): + self.v_bsl[i] = h.GetBinContent(i+1) f.Close() - - - def CorrectBaseline( self ): - """ - apply baseline correction - """ - for pixel in range( self.NPIX ): - self.acalData[pixel,:] -= self.v_bsl[pixel] - - - def info( self ): - """ - print information - """ - print 'data file: ', dfname - print 'calib file: ', calfname + + def baseline_correct(self): + """ subtract baseline from the data + + """ + + for pixel in range(self.npix): + self.acal_data[pixel,:] -= self.v_bsl[pixel] + + def info(self): + """ print run information + + """ + + print 'data file: ', data_file_name + print 'calib file: ', calib_file_name print 'calibration file' - print 'N BaselineMean: ', self.Nblm - print 'N GainMean: ', self.Ngm + print 'N baseline_mean: ', self.Nblm + print 'N gain mean: ', self.Ngm print 'N TriggeroffsetMean: ', self.Ntom -# -------------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- class fnames( object ): """ organize file names of a FACT data run @@ -256,5 +232,5 @@ """ - def __init__( self, specifier = ['012', '023', '2011', '11', '24'], + def __init__(self, specifier = ['012', '023', '2011', '11', '24'], rpath = '/scratch_nfs/res/bsl/', zipped = True): @@ -265,5 +241,7 @@ rpath : directory path for the results; YYYYMMDD will be appended to rpath zipped : use zipped (True) or unzipped (Data) - """ + + """ + self.specifier = specifier self.rpath = rpath @@ -271,5 +249,5 @@ self.make( self.specifier, self.rpath, self.zipped ) - # end of def __init__ + def make( self, specifier, rpath, zipped ): @@ -312,6 +290,4 @@ self.results = self.names['results'] - # end of make - def info( self ): """ print complete filenames @@ -323,61 +299,22 @@ print 'drs-cal: ', self.names['drscal'] print 'results: ', self.names['results'] - # end of def info # end of class definition: fnames( object ) - - -class histogramList( object ): - - def __init__( self, name ): - """ set the name and create empty lists """ - self.name = name # name of the list - self.list = [] # list of the histograms - self.dict = {} # dictionary of histograms - self.hList = TObjArray() # list a la ROOT of the histograms - - def add( self, tag, h ): - self.list.append( h ) - self.dict[tag] = h - self.hList.Add( h ) - - -class pixelHisto1d ( object ): - - def __init__( self, name, title, Nbin, first, last, xtitle, ytitle, NPIX ): - """ - book one dimensional histograms for each pixel - """ - self.name = name - - self.list = [ x for x in range( NPIX ) ] - self.hList = TObjArray() - - for pixel in range( NPIX ): - - hname = name + ' ' + str( pixel ) - htitle = title + ' ' + str( pixel ) - self.list[pixel] = TH1F( hname, htitle, Nbin, first, last ) - - self.list[pixel].GetXaxis().SetTitle( xtitle ) - self.list[pixel].GetYaxis().SetTitle( ytitle ) - self.hList.Add( self.list[pixel] ) - -# simple test method +def _test_iter(): + """ test for function __iter__ """ + +# data_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_111.fits.gz' +# calib_file_name = data_file_name = + data_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_028.fits.gz' + calib_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_022.drs.fits.gz' + run = RawData( data_file_name, calib_file_name ) + + for data, scell, tt in run: + print 'data[0,0] = ', data[0,0], "start_cell[0] =", scell[0], "trigger type = ", tt + + if __name__ == '__main__': - """ - create an instance - """ - dfname = '/data03/fact-construction/raw/2011/11/24/20111124_121.fits' - calfname = '/data03/fact-construction/raw/2011/11/24/20111124_111.drs.fits' - rd = rawdata( dfname, calfname ) - rd.info() - rd.next() - -# for i in range(10): -# df.GetNextRow() - -# print 'evNum: ', evNum.value -# print 'startCells[0:9]: ', startCells[0:9] -# print 'evData[0:9]: ', evData[0:9] + """ tests """ + + _test_iter()