Changeset 17972 for fact/tools/pyscripts/new_pyfact/pyfact.py

Timestamp:

09/24/14 12:45:47 (10 years ago)

Author:

dneise

Message:

added convenience class AuxFile

File:

• fact/tools/pyscripts/new_pyfact/pyfact.py (modified) (12 diffs)

fact/tools/pyscripts/new_pyfact/pyfact.py

@@ -8 +8 @@
 import numpy as np
 import ROOT
-
+import pylab
+
+from scipy.signal import firwin, iirdesign, lfilter
 
 ########## BUILDING OF THE SHARED OBJECT FILES ###############################
@@ -39 +41 @@
     To just build the shared object libs call:
     python pyfact.py build
-    
+
     To build (if necessary) and open an example file
     python -i pyfact.py /path/to/data_file.zfits /path/to/calib_file.drs.fits.gz
-    
+
     Any of the 3 file 'types': fits.gz    zfits    fits
     should be supported.
-    
+
     To clean all of automatically built files, do something like:
     rm *.so *.d AutoDict_* extern_Mars_mcore/*.so extern_Mars_mcore/*.d
@@ -68 +70 @@
 class Fits( object ):
     """ General FITS file access
-    
-    Wrapper for factfits class from Mars/mcore/factfits.h 
+
+    Wrapper for factfits class from Mars/mcore/factfits.h
     (factfits might not be suited to read any FITS file out there, but certainly
     all FACT FITS files can be read using this class)
@@ -82 +84 @@
         self._make_header()
         self._setup_columns()
-    
+
     def _make_header(self):
         """
@@ -97 +99 @@
                 raise IOError("Error: entry type unknown.\n Is %s, but should be one of: [I,F,T,B]" % (entry.type) )
             self.header_comments[key] = entry.comment
-    
+
     def _setup_columns(self):
         """
@@ -113 +115 @@
 
     def next(self, row=None):
-        """ 
+        """
         """
         if row is None:
@@ -124 +126 @@
         return self
 
-    
-    
+class AuxFile( Fits ):
+    """ easy(?) access to FACT aux files
+    """
+    __module__ = 'pyfact'
+    def __init__(self, path, verbose=False):
+        self._verbose = verbose
+        super(AuxFile, self).__init__(path)
+
+    def _setup_columns(self):
+        col_type_to_np_type_map = { 'L' : 'b1',  'A' : 'a1', 'B' : 'i1',
+            'I' : 'i2', 'J' : 'i4', 'K' : 'i8', 'E' : 'f4', 'D' : 'f8'}
+        self._cols = {}
+        self.cols = {}
+        N = self.header['NAXIS2']
+        for key,col in self.f.GetColumns():
+            self._cols[key] = np.zeros(col.num, col_type_to_np_type_map[col.type])
+            self.cols[key] = np.zeros((N, col.num), col_type_to_np_type_map[col.type])
+            if col.num != 0:
+                self.f.SetPtrAddress(key, self._cols[key])
+
+
+        for i,row in enumerate(self):
+            if self._verbose:
+                try:
+                    step = int(self._verbose)
+                except:
+                    step = 10
+                if i % step == 0:
+                    print "reading line", i
+
+            for key in self._cols:
+                self.cols[key][i,:] = self._cols[key]
+
+
+
+
 class RawData( Fits ):
     """ Special raw data FITS file access (with DRS4 calibration)
-        
+
         During iteration the C++ method DrsCalibration::Apply is being called.
     """
@@ -144 +180 @@
         self.drs_calibration = ROOT.DrsCalibration()
         self.drs_calibration.ReadFitsImp( calib_path )
-        
+
         self.drs_calibrate = ROOT.DrsCalibrate()
         self.list_of_previous_start_cells = []
@@ -153 +189 @@
 
     def next(self, row=None):
-        """ 
+        """
         """
         super(RawData, self).next(row)
-        
-        self.drs_calibration.Apply( self.cols['CalibData'], 
-                                    self.cols['Data'], 
-                                    self.cols['StartCellData'], 
+
+        self.drs_calibration.Apply( self.cols['CalibData'],
+                                    self.cols['Data'],
+                                    self.cols['StartCellData'],
                                     self.header['NROI'])
-        
+
         for previous_start_cells in self.list_of_previous_start_cells:
             self.drs_calibrate.CorrectStep(
-                                self.cols['CalibData'], 
-                                self.header['NPIX'], 
-                                self.header['NROI'], 
-                                previous_start_cells, 
-                                self.cols['StartCellData'], 
+                                self.cols['CalibData'],
+                                self.header['NPIX'],
+                                self.header['NROI'],
+                                previous_start_cells,
+                                self.cols['StartCellData'],
                                 self.header['NROI']+10)
             self.drs_calibrate.CorrectStep(
-                                self.cols['CalibData'], 
-                                self.header['NPIX'], 
-                                self.header['NROI'], 
-                                previous_start_cells, 
-                                self.cols['StartCellData'], 
+                                self.cols['CalibData'],
+                                self.header['NPIX'],
+                                self.header['NROI'],
+                                previous_start_cells,
+                                self.cols['StartCellData'],
                                 3)
         self.list_of_previous_start_cells.append(self.cols['StartCellData'])
@@ -192 +228 @@
     print "Example for calibrated raw-file"
     f = RawData(sys.argv[1], sys.argv[2])
-    
+
     print "number of events:", f.header['NAXIS2']
     print "date of observation:", f.header['DATE']
     print "The files has these cols:", f.cols.keys()
-    
+
     for counter,row in enumerate(f):
         print "Event Id:", row.cols['EventNum']
@@ -212 +248 @@
     f.next(3)
     print "Event Id:", f.cols['EventNum']
-    
+
     import matplotlib.pyplot as plt
     plt.ion()
+    """
     for i in range(f.header['NPIX']):
         plt.cla()
@@ -223 +260 @@
         if 'q' in answer.lower():
             break
-
+    """
+
+
+
+    from scipy import signal
+    from scipy.signal import freqs, iirfilter
+    #b, a = iirfilter(4, [1, 10], 1, 60, analog=True, ftype='cheby1')
+    #b, a = signal.iirfilter(17, [50, 200], rs=60, btype='band',
+    #                                analog=True, ftype='cheby2')
+    #w, h = signal.freqs(b, a, 1000)
+    #w, h = freqs(b, a, worN=np.logspace(-1, 2, 1000))
+
+    b, a = iirfilter(N=4, # the order of the filter
+            Wn = [0.15, 0.25], # critical frequencies
+            rp = 1,       # maximum ripple in passband
+            rs = 60,      # minimum attenuation in stopband
+            btype='bandstop', # {'bandpass', 'lowpass', 'highpass', 'bandstop'}
+            #analog=True,
+            ftype='cheby1')
+    #w, h = freqs(b, a, worN=np.logspace(-1, 5, 1000))
+    w, h = freqs(b, a)
+    plt.figure()
+    plt.semilogx(w, 20 * np.log10(abs(h)))
+    plt.xlabel('Frequency')
+    plt.ylabel('Amplitude response [dB]')
+    plt.grid()
+    plt.show()
+
+
+
+
+
+    fig1, (ax11,ax12) = plt.subplots(nrows=2)
+    fig2, (ax21,ax22) = plt.subplots(nrows=2)
+
+    flt = firwin(13, [0.15,0.25])
+    t = np.arange(f.cols['CalibData2D'].shape[1])/2.*1e-9
+    fig3, (ax31, ax32) = plt.subplots(nrows=2)
+    ax31.plot(flt)
+    #w, h = freqs(flt, [1.0])
+    w, h = freqs(flt, [1.0], worN=np.logspace(-1, 5, 1000))
+    ax32.semilogx(w, 20 * np.log10(abs(h)))
+    ax32.set_xlabel('Frequency')
+    ax32.set_ylabel('Amplitude response [dB]')
+    ax32.grid()
+
+
+
+    for i in range(f.header['NPIX']):
+        data = f.cols['CalibData2D'][i]
+        fil_data = lfilter(flt, [1.0], data)
+
+        ax11.cla()
+        ax12.cla()
+        ax11.set_title("Event %d"%(f.cols['EventNum']))
+        ax11.plot(t, data, '.:', label='pixel %d'%i)
+        ax11.legend()
+        ax11.set_ylim(-20,40)
+        ax12.plot(t[:fil_data.size], fil_data, '.:', label='pixel %d'%i)
+        ax12.legend()
+        ax12.set_ylim(-20,40)
+
+        """
+        Pxx, freqs, t, plot = pylab.specgram( f.cols['CalibData2D'][i],
+                NFFT=32, Fs=2000000000,
+                detrend=pylab.detrend_none,
+                window=pylab.window_hanning,
+                noverlap=16)
+        """
+        ax21.cla()
+        ax22.cla()
+        dat, freqs, bins, im = ax21.specgram(data,
+                NFFT=64, Fs=2000000000,
+                detrend=pylab.detrend_none,
+                window=pylab.window_hanning,
+                noverlap=60,
+                interpolation='nearest')
+        ax21.axis('tight')
+        ax21.set_ylabel("frequency [Hz]")
+        ax21.set_xlabel("time [s]")
+
+        dat2, freqs2, bins2, im2 = ax22.specgram(data,
+                NFFT=64, Fs=2000000000,
+                detrend=pylab.detrend_none,
+                window=pylab.window_hanning,
+                noverlap=60,
+                interpolation='nearest')
+        ax22.axis('tight')
+        ax22.set_ylabel("frequency [Hz]")
+        ax22.set_xlabel("time [s]")
+
+        plt.show()
+
+        answer = raw_input('anykey for next pixel; "q" to quit. :')
+        if 'q' in answer.lower():
+            break