Index: fact/tools/pyscripts/new_pyfact/pyfact.py =================================================================== --- fact/tools/pyscripts/new_pyfact/pyfact.py (revision 17795) +++ fact/tools/pyscripts/new_pyfact/pyfact.py (revision 17795) @@ -0,0 +1,218 @@ +#!/usr/bin/python -tt +# +# Werner Lustermann, Dominik Neise +# ETH Zurich, TU Dortmund +# +import os +import ctypes +from ctypes import * +import numpy as np +import pprint # for SlowData +from scipy import signal + +import ROOT +hostname = ROOT.gSystem.HostName() +if 'isdc' in hostname: + ROOT.gSystem.Load("/usr/lib64/libz.so") +elif ('neiseLenovo' in hostname or 'factcontrol' in hostname): + ROOT.gSystem.Load("/usr/lib/libz.so") +elif ("max-K50AB" in hostname or "watz" in hostname): + ROOT.gSystem.Load("/usr/lib/x86_64-linux-gnu/libz.so") +elif ("grolsch" in hostname): + ROOT.gSystem.Load("/usr/lib/i386-linux-gnu/libz.so") +else: + print "Error,Warning,Whatever libz stuff makes me crazy." + +root_make_string = ROOT.gSystem.GetMakeSharedLib() +if not "-std=c++0x" in root_make_string: + make_string = root_make_string.replace('$Opt', '$Opt -std=c++0x -D HAVE_ZLIB') +ROOT.gSystem.SetMakeSharedLib(make_string) +ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/izstream.h+O") +ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/fits.h+O") +ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/zfits.h+O") +ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/factfits.h+O") +if not "-std=c++0x" in root_make_string: + make_string = root_make_string.replace('$Opt', "$Opt -std=c++0x -D HAVE_ZLIB -D'PACKAGE_NAME=\"PACKAGE_NAME\"' -D'PACKAGE_VERSION=\"PACKAGE_VERSION\"' -D'REVISION=\"REVISION\"' ") +ROOT.gSystem.SetMakeSharedLib(make_string) +ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/DrsCalib.h+O") + +ROOT.gInterpreter.GenerateDictionary("map","map;string;extern_Mars_mcore/fits.h") +ROOT.gInterpreter.GenerateDictionary("pair","map;string;extern_Mars_mcore/fits.h") +ROOT.gInterpreter.GenerateDictionary("map","map;string;extern_Mars_mcore/fits.h") +ROOT.gInterpreter.GenerateDictionary("pair","map;string;extern_Mars_mcore/fits.h") +ROOT.gInterpreter.GenerateDictionary("vector","vector;extern_Mars_mcore/DrsCalib.h") + +ROOT.gSystem.Load('extern_Mars_mcore/fits_h.so') +ROOT.gSystem.Load('extern_Mars_mcore/izstream_h.so') +ROOT.gSystem.Load('extern_Mars_mcore/zfits_h.so') +ROOT.gSystem.Load('extern_Mars_mcore/factfits_h.so') +ROOT.gSystem.Load('extern_Mars_mcore/DrsCalib_h.so') +from ROOT import * + +class Fits( object ): + """ General FITS file access + + 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) + """ + __module__ = 'pyfact' + def __init__(self, path): + """ + """ + if not os.path.exists(path): + raise IOError(path+' was not found') + self.f = factfits(path) + self._make_header() + self._setup_columns() + + def _make_header(self): + """ + """ + str_to_bool = { 'T':True, 'F':False} + type_conversion = { 'I' : int, 'F' : float, 'T' : str, 'B' : str_to_bool.__getitem__} + + self.header = {} + self.header_comments = {} + for key,entry in self.f.GetKeys(): + try: + self.header[key] = type_conversion[entry.type](entry.value) + except KeyError: + 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): + """ + """ + col_type_to_np_type_map = { 'L' : 'b1', 'A' : 'a1', 'B' : 'i1', + 'I' : 'i2', 'J' : 'i4', 'K' : 'i8', 'E' : 'f4', 'D' : 'f8'} + self.cols = {} + for key,col in self.f.GetColumns(): + self.cols[key] = np.zeros(col.num, col_type_to_np_type_map[col.type]) + if col.num != 0: + self.f.SetPtrAddress(key, self.cols[key]) + + def __iter__(self): + return self + + def next(self, row=None): + """ + """ + if row is None: + if self.f.GetNextRow() == False: + raise StopIteration + else: + row = int(row) + if self.f.GetRow(row) == False: + raise StopIteration + return self + + + +class RawData( Fits ): + """ Special raw data FITS file access (with DRS4 calibration) + + During iteration the C++ method DrsCalibration::Apply is being called. + """ + __module__='pyfact' + def __init__(self, data_path, calib_path): + """ -constructor- + *data_path* : fits or fits.gz file of the data including the path + *calib_path* : fits or fits.gz file containing DRS calibration data + """ + super(RawData, self).__init__(data_path) + self.cols['CalibData'] = np.zeros( self.cols['Data'].shape, np.float32) + self.cols['CalibData2D'] = self.cols['CalibData'].reshape( self.header['NPIX'], -1) + if not self.cols['CalibData2D'].base is self.cols['CalibData']: + print "Error seomthing went wrong!" + self.drs_calibration = DrsCalibration() + self.drs_calibration.ReadFitsImp( calib_path ) + + self.drs_calibrate = DrsCalibrate() + self.list_of_previous_start_cells = [] + + def __iter__(self): + """ iterator """ + return self + + def next(self, row=None): + """ + """ + super(RawData, self).next(row) + + 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.header['NROI']+10) + ) + self.drs_calibrate.CorrectStep( + self.cols['CalibData'], + self.header['NPIX'], + self.header['NROI'], + previous_start_cells, + self.cols['StartCellData'], + 3) + ) + self.fPrevStart.append(self.cols['StartCellData']) + if len(self.fPrevStart) > 5: + self.fPrevStart.pop(0) + + for ch in range(self.header['NPIX']): + self.drs_calibrate.RemoveSpikes3(self.cols['CalibData2D'][ch], self.header['NROI']) + + return self + + + +if __name__ == '__main__': + """ tests """ + import sys + if len(sys.argv) == 2: + print "Example for aux-file or uncalibrated raw-file" + f = Fits(sys.argv[1]) + elif len(sys.argv) == 3: + print "Example for calibrated raw-file" + f = RawData(sys.argv[1], sys.argv[2]) + else: + f = None + + 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'] + print "shape of column 'StartCellData'", row.cols['StartCellData'].shape + print "dtype of column 'Data'", row.cols['StartCellData'].dtype + if counter > 3: + break + # get next row + f.next() + print "Event Id:", f.cols['EventNum'] + # get another row + f.next(10) + print "Event Id:", f.cols['EventNum'] + # Go back again + f.next(3) + print "Event Id:", f.cols['EventNum'] + + if len(sys.argv) == 3: + import matplotlib.pyplot as plt + plt.ion() + for i in range(f.header['NPIX']): + plt.cla() + plt.plot(f.cols['CalibData2D'][i], '.:', label='pixel %d'%i) + plt.legend() + answer = raw_input('anykey for next pixel; "q" to quit. :') + if 'q' in answer.lower(): + break +