Index: fact/tools/pyscripts/new_pyfact/pyfact.py =================================================================== --- fact/tools/pyscripts/new_pyfact/pyfact.py (revision 17915) +++ fact/tools/pyscripts/new_pyfact/pyfact.py (revision 17972) @@ -8,5 +8,7 @@ import numpy as np import ROOT - +import pylab + +from scipy.signal import firwin, iirdesign, lfilter ########## BUILDING OF THE SHARED OBJECT FILES ############################### @@ -39,11 +41,11 @@ 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,6 +70,6 @@ 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,5 +84,5 @@ self._make_header() self._setup_columns() - + def _make_header(self): """ @@ -97,5 +99,5 @@ 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,5 +115,5 @@ def next(self, row=None): - """ + """ """ if row is None: @@ -124,9 +126,43 @@ 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,5 +180,5 @@ self.drs_calibration = ROOT.DrsCalibration() self.drs_calibration.ReadFitsImp( calib_path ) - + self.drs_calibrate = ROOT.DrsCalibrate() self.list_of_previous_start_cells = [] @@ -153,27 +189,27 @@ 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,9 +228,9 @@ 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,7 +248,8 @@ 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,3 +260,98 @@ 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