#!/usr/bin/python # # Werner Lustermann # ETH Zurich # from ctypes import * import numpy as np from scipy import signal # get the ROOT stuff + my shared libs from ROOT import gSystem # fitslib.so is made from fits.h and is used to access the data gSystem.Load('~/py/fitslib.so') from ROOT import * class rawdata( object ): """ raw data access and calibration - open raw data file and drs calibration file - performs amplitude calibration - performs baseline substraction if wanted - provides all data in an array: row = number of pixel col = length of region of interest """ # constructor of the classe def __init__( self, dfname, calfname, bslfname='' ): """ open data file and calibration data file get basic information about the data in dfname 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 # baseline correction: True / False if len( bslfname ) == 0: self.correct_baseline = False else: self.correct_baseline = True # access data file try: df = fits( self.dfname ) except IOError: print 'problem accessing data file: ', dfname raise # stop ! no data self.df = df # 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 # 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 ) # 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 try: calf = fits( self.calfname ) except IOError: print 'problem accessing calibration file: ', calfname 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 ): # 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] 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' ): """ 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 ) 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 ) 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 print 'calibration file' print 'N BaselineMean: ', self.Nblm print 'N GainMean: ', self.Ngm print 'N TriggeroffsetMean: ', self.Ntom # -------------------------------------------------------------------------------- class fnames( object ): """ organize file names of a FACT data run """ def __init__( self, specifier = ['012', '023', '2011', '11', '24'], rpath = '/scratch_nfs/bsl/', zipped = True): """ specifier : list of strings defined as: [ 'DRS calibration file', 'Data file', 'YYYY', 'MM', 'DD'] 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 self.zipped = zipped self.make( self.specifier, self.rpath, self.zipped ) # end of def __init__ def make( self, specifier, rpath, zipped ): """ create (make) the filenames names : dictionary of filenames, tags { 'data', 'drscal', 'results' } data : name of the data file drscal : name of the drs calibration file results : radikal of file name(s) for results (to be completed by suffixes) """ self.specifier = specifier if zipped: dpath = '/data00/fact-construction/raw/' ext = '.fits.gz' else: dpath = '/data03/fact-construction/raw/' ext = '.fits' year = specifier[2] month = specifier[3] day = specifier[4] yyyymmdd = year + month + day dfile = specifier[1] cfile = specifier[0] rpath = rpath + yyyymmdd + '/' self.rpath = rpath self.names = {} tmp = dpath + year + '/' + month + '/' + day + '/' + yyyymmdd + '_' self.names['data'] = tmp + dfile + ext self.names['drscal'] = tmp + cfile + '.drs' + ext self.names['results'] = rpath + yyyymmdd + '_' + dfile + '_' + cfile self.data = self.names['data'] self.drscal = self.names['drscal'] self.results = self.names['results'] # end of make def info( self ): """ print complete filenames """ print 'file names:' print 'data: ', self.names['data'] 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 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]