Changeset 13124 for fact/tools


Ignore:
Timestamp:
03/15/12 22:38:09 (13 years ago)
Author:
lusterma
Message:
pyfact made identical to pyfact_rename.py
File:
1 edited

Legend:

Unmodified
Added
Removed

  • fact/tools/pyscripts/pyfact/pyfact.py

    r12947 r13124  
    44# ETH Zurich
    55#
    6 from   ctypes import *
     6from ctypes import *
    77import numpy as np
    88from scipy import signal
    … …  
    1414from ROOT import *
    1515
    16 class rawdata( object ):
    17     """
    18     raw data access and calibration
     16
     17class RawData( object ):
     18    """ raw data access and calibration
     19   
    1920    - open raw data file and drs calibration file
    2021    - performs amplitude calibration
    … …  
    2324      row = number of pixel
    2425      col = length of region of interest
     26     
    2527    """
    26     # constructor of the classe
    27     def __init__( self, dfname,  calfname, bslfname='' ):
    28         """
     28
     29    def __init__(self, data_file_name, calib_file_name,
     30                 user_action_calib=lambda acal_data, data, blm, tom, gm, scells, nroi: None,
     31                 baseline_file_name=''):
     32        """ initialize object
     33
    2934        open data file and calibration data file
    30         get basic information about the data in dfname
     35        get basic information about the data in data_file_name
    3136        allocate buffers for data access
    3237
    33         dfname   : fits or fits.gz file containing the data including the path
    34         calfname : fits or fits.gz file containing DRS calibration data
    35         bslfname : npy file containing the baseline values
    36         """
    37         self.dfname   = dfname
    38         self.calfname = calfname
    39         self.bslfname = bslfname
     38        data_file_name   : fits or fits.gz file of the data including the path
     39        calib_file_name : fits or fits.gz file containing DRS calibration data
     40        baseline_file_name : npy file containing the baseline values
     41       
     42        """
     43
     44        self.data_file_name = data_file_name
     45        self.calib_file_name = calib_file_name
     46        self.baseline_file_name = baseline_file_name
     47           
     48        self.user_action_calib = user_action_calib
    4049       
    4150        # baseline correction: True / False
    42         if len( bslfname ) == 0:
     51        if len(baseline_file_name) == 0:
    4352            self.correct_baseline = False
    4453        else:
    … …  
    4756        # access data file
    4857        try:
    49             df = fits( self.dfname )
     58            data_file = fits(self.data_file_name)
    5059        except IOError:
    51             print 'problem accessing data file: ', dfname
    52             raise # stop ! no data
    53         self.df = df
     60            print 'problem accessing data file: ', data_file_name
     61            raise  # stop ! no data
     62        #: data file (fits object)
     63        self.data_file = data_file
    5464       
    5565        # get basic information about the data file
    56         self.NROI    = df.GetUInt( 'NROI' ) # region of interest (length of DRS pipeline read out)
    57         self.NPIX    = df.GetUInt( 'NPIX' ) # number of pixels (should be 1440)
    58         self.NEvents = df.GetNumRows()      # find number of events
     66        #: region of interest (number of DRS slices read)
     67        self.nroi    = data_file.GetUInt('NROI')
     68        #: number of pixels (should be 1440)
     69        self.npix    = data_file.GetUInt('NPIX')
     70        #: number of events in the data run
     71        self.nevents = data_file.GetNumRows()
     72       
    5973        # allocate the data memories
    60         self.evNum = c_ulong()
    61         self.trigType = c_ushort()
    62         self.Data  = np.zeros( self.NPIX * self.NROI, np.int16 )
    63         self.startCells = np.zeros( self.NPIX, np.int16 )
     74        self.event_id = c_ulong()
     75        self.trigger_type = c_ushort()
     76        #: 1D array with raw data
     77        self.data  = np.zeros( self.npix * self.nroi, np.int16 )
     78        #: slice where drs readout started
     79        self.start_cells = np.zeros( self.npix, np.int16 )
     80        #: time when the FAD was triggered, in some strange units...
     81        self.board_times = np.zeros( 40, np.int32 )
     82
    6483        # set the pointers to the data++
    65         df.SetPtrAddress( 'EventNum', self.evNum )
    66         df.SetPtrAddress( 'TriggerType', self.trigType )
    67         df.SetPtrAddress( 'StartCellData', self.startCells ) # DRS readout start cell
    68         df.SetPtrAddress( 'Data', self.Data ) # this is what you would expect
    69         # df.GetNextRow() # access the first event
    70        
    71         # access calibration file
     84        data_file.SetPtrAddress('EventNum', self.event_id)
     85        data_file.SetPtrAddress('TriggerType', self.trigger_type)
     86        data_file.SetPtrAddress('StartCellData', self.start_cells)
     87        data_file.SetPtrAddress('Data', self.data)
     88        data_file.SetPtrAddress('BoardTime', self.board_times)
     89               
     90        # open the calibration file
    7291        try:
    73             calf = fits( self.calfname )
     92            calib_file = fits(self.calib_file_name)
    7493        except IOError:
    75             print 'problem accessing calibration file: ', calfname
     94            print 'problem accessing calibration file: ', calib_file_name
    7695            raise
    77         self.calf = calf
    78         #
    79         BaselineMean      = calf.GetN('BaselineMean')
    80         GainMean          = calf.GetN('GainMean')
    81         TriggerOffsetMean = calf.GetN('TriggerOffsetMean')
    82 
    83         self.blm = np.zeros( BaselineMean, np.float32 )
    84         self.gm  = np.zeros( GainMean, np.float32 )
    85         self.tom = np.zeros( TriggerOffsetMean, np.float32 )
    86 
    87         self.Nblm = BaselineMean / self.NPIX
    88         self.Ngm  = GainMean / self.NPIX
    89         self.Ntom  = TriggerOffsetMean / self.NPIX
    90 
    91         calf.SetPtrAddress( 'BaselineMean', self.blm )
    92         calf.SetPtrAddress( 'GainMean', self.gm )
    93         calf.SetPtrAddress( 'TriggerOffsetMean', self.tom )
    94         calf.GetRow(0)
    95 
    96         self.v_bsl = np.zeros( self.NPIX ) # array with baseline values (all ZERO)
    97         self.smoothData = None
    98         self.maxPos = None
    99         self.maxAmp = None
    100 
    101 
    102     def next( self ):
    103         """
    104         load the next event from disk and calibrate it
    105         """
    106         self.df.GetNextRow()
    107         self.calibrate_drsAmplitude()
    108 
    109        
    110     def calibrate_drsAmplitude( self ):
    111         """
    112         perform amplitude calibration for the event
    113         """
    114         tomV = 2000./4096.
    115         acalData = self.Data * tomV # convert into mV
    116 
    117         # reshape arrays: row = pixel, col = drs_slice
    118         acalData = np.reshape( acalData, (self.NPIX, self.NROI) )
    119         blm = np.reshape( self.blm, (self.NPIX, 1024) )
    120         tom = np.reshape( self.tom, (self.NPIX, 1024) )
    121         gm  = np.reshape( self.gm,  (self.NPIX, 1024) )
    122        
    123         # print 'acal Data ', acalData.shape
    124         # print 'blm shape ', blm.shape
    125         # print 'gm shape  ', gm.shape
    126        
    127         for pixel in range( self.NPIX ):
     96        #: drs calibration file
     97        self.calib_file = calib_file
     98       
     99        baseline_mean = calib_file.GetN('BaselineMean')
     100        gain_mean = calib_file.GetN('GainMean')
     101        trigger_offset_mean = calib_file.GetN('TriggerOffsetMean')
     102
     103        self.blm = np.zeros(baseline_mean, np.float32)
     104        self.gm  = np.zeros(gain_mean, np.float32)
     105        self.tom = np.zeros(trigger_offset_mean, np.float32)
     106
     107        self.Nblm = baseline_mean / self.npix
     108        self.Ngm  = gain_mean / self.npix
     109        self.Ntom  = trigger_offset_mean / self.npix
     110
     111        calib_file.SetPtrAddress('BaselineMean', self.blm)
     112        calib_file.SetPtrAddress('GainMean', self.gm)
     113        calib_file.SetPtrAddress('TriggerOffsetMean', self.tom)
     114        calib_file.GetRow(0)
     115
     116        self.v_bsl = np.zeros(self.npix)  # array of baseline values (all ZERO)
     117        self.data_saverage_out = None
     118        self.pulse_time_of_maximum = None
     119        self.pulse_amplitude = None
     120
     121    def __iter__(self):
     122        """ iterator """
     123        return self
     124       
     125    def __add__(self, jump_over):
     126        self.data_file.GetRow(jump_over)
     127        return self
     128       
     129    def next(self):
     130        """ used by __iter__ """
     131
     132        if self.data_file.GetNextRow() == False:
     133            raise StopIteration
     134        else:
     135            self.calibrate_drs_amplitude()
     136
     137        #print 'nevents = ', self.nevents, 'event_id = ', self.event_id.value
     138       
     139        return self.acal_data, self.start_cells, self.trigger_type.value
     140
     141    def next_event(self):
     142        """ load the next event from disk and calibrate it
     143       
     144        """
     145
     146        self.data_file.GetNextRow()
     147        self.calibrate_drs_amplitude()
     148
     149    def calibrate_drs_amplitude(self):
     150        """ perform the drs amplitude calibration of the event data
     151       
     152        """
     153
     154        to_mV = 2000./4096.
     155        #: 2D array with amplitude calibrated dat in mV
     156        acal_data = self.data * to_mV  # convert ADC counts to mV
     157
     158        # make 2D arrays: row = pixel, col = drs_slice
     159        acal_data = np.reshape(acal_data, (self.npix, self.nroi) )
     160        blm = np.reshape(self.blm, (self.npix, self.Nblm) )
     161        tom = np.reshape(self.tom, (self.npix, self.Ntom) )
     162        gm  = np.reshape(self.gm,  (self.npix, self.Ngm) )
     163       
     164        for pixel in range( self.npix ):
    128165            # rotate the pixel baseline mean to the Data startCell
    129             blm_pixel = np.roll( blm[pixel,:], -self.startCells[pixel] )
    130             acalData[pixel,:] -= blm_pixel[0:self.NROI]
    131             acalData[pixel,:] -= tom[pixel, 0:self.NROI]
    132             acalData[pixel,:] /= gm[pixel,  0:self.NROI]
     166            blm_pixel = np.roll( blm[pixel,:], -self.start_cells[pixel] )
     167            tom_pixel = np.roll( tom[pixel,:], -self.start_cells[pixel] )
     168            gm_pixel = np.roll( gm[pixel,:], -self.start_cells[pixel] )
     169            acal_data[pixel,:] -= blm_pixel[0:self.nroi]
     170            acal_data[pixel,:] -= tom_pixel[0:self.nroi]
     171            acal_data[pixel,:] /= gm_pixel[0:self.nroi]
    133172           
    134         self.acalData = acalData * 1907.35
    135    
    136         # print 'acalData ', self.acalData[0:2,0:20]
    137        
    138     def filterSlidingAverage( self , windowSize = 4):
    139         """ sliding average filter
    140         using:
    141             self.acalData
    142         filling array:
    143             self.smoothData
    144         """
    145         #scipy.signal.lfilter(b, a, x, axis=-1, zi=None)
    146         smoothData = self.acalData.copy()
    147         b = np.ones( windowSize )
    148         a = np.zeros( windowSize )
    149         a[0] = len(b)
    150         smoothData[:,:] = signal.lfilter(b, a, smoothData[:,:])
    151 
    152         self.smoothData = smoothData
    153        
    154     def filterCFD( self, length=10, ratio=0.75):
    155         """ constant fraction filter
    156         using:
    157             self.smoothData
    158         filling array:
    159             self.cfdData
    160         """
    161         if self.smoothData == None:
    162             print 'error pyfact.filterCFD was called without prior call to filterSlidingAverage'
    163             print ' variable self.smoothData is needed '
    164             pass
    165         cfdData = self.smoothData.copy()
    166         b = np.zeros( length )
    167         a = np.zeros( length )
    168         b[0] = -1. * ratio
    169         b[length-1] = 1.
    170         a[0] = 1.
    171         cfdData[:,:] = signal.lfilter(b, a, cfdData[:,:])
    172        
    173         self.cfdData = cfdData
    174    
    175     def findPeak (self, min=30, max=250):
    176         """ find maximum in search window
    177         using:
    178             self.smoothData
    179         filling arrays:
    180             self.maxPos
    181             self.maxAmp
    182         """
    183         if self.smoothData == None:
    184             print 'error pyfact.findPeakMax was called without prior call to filterSlidingAverage'
    185             print ' variable self.smoothData is needed '
    186             pass
    187         maxPos = np.argmax( self.smoothData[:,min:max] , 1)
    188         maxAmp = np.max( self.smoothData[:,min:max] , 1)
    189         self.maxPos = maxPos
    190         self.maxAmp = maxAmp
    191 
    192     def sumAroundPeak (self, left=13, right=23):
    193         """ integrate signal in gate around Peak
    194         using:
    195             self.maxPos
    196             self.acalData
    197         filling array:
    198             self.integral
    199         """
    200         if self.maxPos == None:
    201             print 'error pyfact.sumAroundPeak was called without prior call of findPeak'
    202             print ' variable self.maxPos is needed'
    203             pass
    204        
    205         sums = np.empty( self.NPIX )
    206         for pixel in range( self.NPIX ):
    207             min = self.maxPos[pixel]-left
    208             max = self.maxPos[pixel]+right
    209             sums[pixel] = self.acalData[pixel,min:max].sum()
    210        
    211         self.integral = sums
    212        
    213     def ReadBaseline( self, file, bsl_hist = 'bsl_sum/hplt_mean' ):
    214         """
     173        self.acal_data = acal_data * 1907.35
     174       
     175        #print 'blm _pyfact', blm[0,0:20]
     176        #t = np.roll( blm[0,:], -self.start_cells[0] )
     177        #print 'blm _pyfact', t[0:20]
     178        #print 'start_pyfact: ', self.start_cells[0]
     179        #print 'acal _pyfact: ', self.acal_data[0,0:10]
     180        #t = np.roll( gm[0,:], -self.start_cells[0] )
     181        #print 'gm _pyfact: ', t[0:10]
     182        self.user_action_calib( self.acal_data,
     183                np.reshape(self.data, (self.npix, self.nroi) ), blm, tom, gm, self.start_cells, self.nroi)
     184
     185       
     186    def baseline_read_values(self, file, bsl_hist='bsl_sum/hplt_mean'):
     187        """
     188       
    215189        open ROOT file with baseline histogram and read baseline values
    216190        file       name of the root file
    217191        bsl_hist   path to the histogram containing the basline values
    218         """
     192
     193        """
     194
    219195        try:
    220             f = TFile( file )
     196            f = TFile(file)
    221197        except:
    222198            print 'Baseline data file could not be read: ', file
    223199            return
    224200       
    225         h = f.Get( bsl_hist )
    226 
    227         for i in range( self.NPIX ):
    228             self.v_bsl[i] = h.GetBinContent( i+1 )
     201        h = f.Get(bsl_hist)
     202
     203        for i in range(self.npix):
     204            self.v_bsl[i] = h.GetBinContent(i+1)
    229205
    230206        f.Close()
    231 
    232        
    233     def CorrectBaseline( self ):
    234         """
    235         apply baseline correction
    236         """
    237         for pixel in range( self.NPIX ):
    238             self.acalData[pixel,:] -= self.v_bsl[pixel]
    239            
    240        
    241     def info( self ):
    242         """
    243         print information
    244         """
    245         print 'data file:  ', dfname
    246         print 'calib file: ', calfname
     207       
     208    def baseline_correct(self):
     209        """ subtract baseline from the data
     210
     211        """
     212       
     213        for pixel in range(self.npix):
     214            self.acal_data[pixel,:] -= self.v_bsl[pixel]
     215                   
     216    def info(self):
     217        """ print run information
     218       
     219        """
     220       
     221        print 'data file:  ', data_file_name
     222        print 'calib file: ', calib_file_name
    247223        print 'calibration file'
    248         print 'N BaselineMean: ', self.Nblm
    249         print 'N GainMean: ', self.Ngm
     224        print 'N baseline_mean: ', self.Nblm
     225        print 'N gain mean: ', self.Ngm
    250226        print 'N TriggeroffsetMean: ', self.Ntom
    251227       
    252 # --------------------------------------------------------------------------------
     228# -----------------------------------------------------------------------------
    253229class fnames( object ):
    254230    """ organize file names of a FACT data run
    … …  
    256232    """
    257233   
    258     def __init__( self, specifier = ['012', '023', '2011', '11', '24'],
     234    def __init__(self, specifier = ['012', '023', '2011', '11', '24'],
    259235                 rpath = '/scratch_nfs/res/bsl/',
    260236                 zipped = True):
    … …  
    265241        rpath     : directory path for the results; YYYYMMDD will be appended to rpath
    266242        zipped    : use zipped (True) or unzipped (Data)
    267         """
     243
     244        """
     245       
    268246        self.specifier = specifier
    269247        self.rpath     = rpath
    … …  
    271249       
    272250        self.make( self.specifier, self.rpath, self.zipped )
    273     # end of def __init__
     251
    274252
    275253    def make( self, specifier, rpath, zipped ):
    … …  
    312290        self.results = self.names['results']
    313291
    314     # end of make
    315 
    316292    def info( self ):
    317293        """ print complete filenames
    … …  
    323299        print 'drs-cal: ', self.names['drscal']
    324300        print 'results: ', self.names['results']
    325     # end of def info
    326301
    327302# end of class definition: fnames( object )
    328303
    329 
    330 
    331 class histogramList( object ):
    332 
    333     def __init__( self, name ):
    334         """ set the name and create empty lists """
    335         self.name  = name         # name of the list
    336         self.list  = []           # list of the histograms
    337         self.dict  = {}           # dictionary of histograms
    338         self.hList = TObjArray()  # list a la ROOT of the histograms
    339 
    340     def add( self, tag, h ):
    341         self.list.append( h )
    342         self.dict[tag] = h
    343         self.hList.Add( h )
    344 
    345 
    346 class pixelHisto1d ( object ):
    347 
    348     def __init__( self, name, title, Nbin, first, last, xtitle, ytitle, NPIX ):
    349         """
    350         book one dimensional histograms for each pixel
    351         """
    352         self.name = name
    353 
    354         self.list = [ x for x in range( NPIX ) ]
    355         self.hList = TObjArray()
    356 
    357         for pixel in range( NPIX ):
    358 
    359             hname  = name + ' ' + str( pixel )
    360             htitle = title + ' ' + str( pixel )
    361             self.list[pixel] = TH1F( hname, htitle, Nbin, first, last )
    362 
    363             self.list[pixel].GetXaxis().SetTitle( xtitle )
    364             self.list[pixel].GetYaxis().SetTitle( ytitle )
    365             self.hList.Add( self.list[pixel] )
    366 
    367 # simple test method
     304def _test_iter():
     305    """ test for function __iter__ """
     306
     307#    data_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_111.fits.gz'
     308#    calib_file_name =     data_file_name =
     309    data_file_name =  '/home/luster/win7/FACT/data/raw/20120114/20120114_028.fits.gz'
     310    calib_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_022.drs.fits.gz'
     311    run = RawData( data_file_name, calib_file_name )
     312
     313    for data, scell, tt in run:
     314        print 'data[0,0] = ', data[0,0], "start_cell[0] =", scell[0], "trigger type = ", tt
     315
     316
    368317if __name__ == '__main__':
    369     """
    370     create an instance
    371     """
    372     dfname = '/data03/fact-construction/raw/2011/11/24/20111124_121.fits'
    373     calfname = '/data03/fact-construction/raw/2011/11/24/20111124_111.drs.fits'
    374     rd = rawdata( dfname, calfname )
    375     rd.info()
    376     rd.next()
    377    
    378 # for i in range(10):
    379 #    df.GetNextRow()
    380 
    381 #    print 'evNum: ', evNum.value
    382 #    print 'startCells[0:9]: ', startCells[0:9]
    383 #    print 'evData[0:9]: ', evData[0:9]
     318    """ tests  """
     319
     320    _test_iter()
Note: See TracChangeset for help on using the changeset viewer.