Ignore:
Timestamp:
02/03/12 09:59:59 (13 years ago)
Author:
neise
Message:
implemented filters, peak-amplitude and peak-integral computation
File:
1 edited

Legend:

Unmodified
Added
Removed
  • fact/tools/pyscripts/pyfact/pyfact.py

    r12815 r12816  
    66from ctypes import *
    77import numpy as np
     8import scipy.signal as spsi
    89
    910# get the ROOT stuff + my shared libs
     
    5859        # allocate the data memories
    5960        self.evNum = c_ulong()
     61        self.trigType = c_ushort()
    6062        self.Data  = np.zeros( self.NPIX * self.NROI, np.int16 )
    6163        self.startCells = np.zeros( self.NPIX, np.int16 )
    6264        # set the pointers to the data++
    6365        df.SetPtrAddress( 'EventNum', self.evNum )
     66        df.SetPtrAddress( 'TriggerType', self.trigType )
    6467        df.SetPtrAddress( 'StartCellData', self.startCells ) # DRS readout start cell
    6568        df.SetPtrAddress( 'Data', self.Data ) # this is what you would expect
     
    100103        self.df.GetNextRow()
    101104        self.calibrate_drsAmplitude()
     105        self.smoothData = None
     106        self.maxPos = None
     107        self.maxAmp = None
    102108
    103109       
     
    129135   
    130136        # print 'acalData ', self.acalData[0:2,0:20]
    131 
     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[:,:] = spsi.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[:,:] = spsi.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       
    132213    def ReadBaseline( self, file, bsl_hist = 'bsl_sum/hplt_mean' ):
    133214        """
Note: See TracChangeset for help on using the changeset viewer.