| 1 | #!/usr/bin/python -tt
|
|---|
| 2 | #
|
|---|
| 3 | # Werner Lustermann, Dominik Neise
|
|---|
| 4 | # ETH Zurich, TU Dortmund
|
|---|
| 5 | #
|
|---|
| 6 | import os
|
|---|
| 7 | import ctypes
|
|---|
| 8 | from ctypes import *
|
|---|
| 9 | import numpy as np
|
|---|
| 10 | import pprint # for SlowData
|
|---|
| 11 | from scipy import signal
|
|---|
| 12 |
|
|---|
| 13 | import ROOT
|
|---|
| 14 | hostname = ROOT.gSystem.HostName()
|
|---|
| 15 | if 'isdc' in hostname:
|
|---|
| 16 | ROOT.gSystem.Load("/usr/lib64/libz.so")
|
|---|
| 17 | elif ('neiseLenovo' in hostname or 'factcontrol' in hostname):
|
|---|
| 18 | ROOT.gSystem.Load("/usr/lib/libz.so")
|
|---|
| 19 | elif ("max-K50AB" in hostname or "watz" in hostname):
|
|---|
| 20 | ROOT.gSystem.Load("/usr/lib/x86_64-linux-gnu/libz.so")
|
|---|
| 21 | elif ("grolsch" in hostname):
|
|---|
| 22 | ROOT.gSystem.Load("/usr/lib/i386-linux-gnu/libz.so")
|
|---|
| 23 | else:
|
|---|
| 24 | print "Error,Warning,Whatever libz stuff makes me crazy."
|
|---|
| 25 |
|
|---|
| 26 | root_make_string = ROOT.gSystem.GetMakeSharedLib()
|
|---|
| 27 | if not "-std=c++0x" in root_make_string:
|
|---|
| 28 | make_string = root_make_string.replace('$Opt', '$Opt -std=c++0x -D HAVE_ZLIB')
|
|---|
| 29 | ROOT.gSystem.SetMakeSharedLib(make_string)
|
|---|
| 30 | ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/izstream.h+O")
|
|---|
| 31 | ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/fits.h+O")
|
|---|
| 32 | ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/zfits.h+O")
|
|---|
| 33 | ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/factfits.h+O")
|
|---|
| 34 | if not "-std=c++0x" in root_make_string:
|
|---|
| 35 | 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\"' ")
|
|---|
| 36 | ROOT.gSystem.SetMakeSharedLib(make_string)
|
|---|
| 37 | ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/DrsCalib.h+O")
|
|---|
| 38 |
|
|---|
| 39 | ROOT.gInterpreter.GenerateDictionary("map<string,fits::Entry>","map;string;extern_Mars_mcore/fits.h")
|
|---|
| 40 | ROOT.gInterpreter.GenerateDictionary("pair<string,fits::Entry>","map;string;extern_Mars_mcore/fits.h")
|
|---|
| 41 | ROOT.gInterpreter.GenerateDictionary("map<string,fits::Table::Column>","map;string;extern_Mars_mcore/fits.h")
|
|---|
| 42 | ROOT.gInterpreter.GenerateDictionary("pair<string,fits::Table::Column>","map;string;extern_Mars_mcore/fits.h")
|
|---|
| 43 | ROOT.gInterpreter.GenerateDictionary("vector<DrsCalibrate::Step>","vector;extern_Mars_mcore/DrsCalib.h")
|
|---|
| 44 |
|
|---|
| 45 | ROOT.gSystem.Load('extern_Mars_mcore/fits_h.so')
|
|---|
| 46 | ROOT.gSystem.Load('extern_Mars_mcore/izstream_h.so')
|
|---|
| 47 | ROOT.gSystem.Load('extern_Mars_mcore/zfits_h.so')
|
|---|
| 48 | ROOT.gSystem.Load('extern_Mars_mcore/factfits_h.so')
|
|---|
| 49 | ROOT.gSystem.Load('extern_Mars_mcore/DrsCalib_h.so')
|
|---|
| 50 | from ROOT import *
|
|---|
| 51 |
|
|---|
| 52 | class Fits( object ):
|
|---|
| 53 | """ General FITS file access
|
|---|
| 54 |
|
|---|
| 55 | Wrapper for factfits class from Mars/mcore/factfits.h
|
|---|
| 56 | (factfits might not be suited to read any FITS file out there, but certainly
|
|---|
| 57 | all FACT FITS files can be read using this class)
|
|---|
| 58 | """
|
|---|
| 59 | __module__ = 'pyfact'
|
|---|
| 60 | def __init__(self, path):
|
|---|
| 61 | """
|
|---|
| 62 | """
|
|---|
| 63 | if not os.path.exists(path):
|
|---|
| 64 | raise IOError(path+' was not found')
|
|---|
| 65 | self.f = factfits(path)
|
|---|
| 66 | self._make_header()
|
|---|
| 67 | self._setup_columns()
|
|---|
| 68 |
|
|---|
| 69 | def _make_header(self):
|
|---|
| 70 | """
|
|---|
| 71 | """
|
|---|
| 72 | str_to_bool = { 'T':True, 'F':False}
|
|---|
| 73 | type_conversion = { 'I' : int, 'F' : float, 'T' : str, 'B' : str_to_bool.__getitem__}
|
|---|
| 74 |
|
|---|
| 75 | self.header = {}
|
|---|
| 76 | self.header_comments = {}
|
|---|
| 77 | for key,entry in self.f.GetKeys():
|
|---|
| 78 | try:
|
|---|
| 79 | self.header[key] = type_conversion[entry.type](entry.value)
|
|---|
| 80 | except KeyError:
|
|---|
| 81 | raise IOError("Error: entry type unknown.\n Is %s, but should be one of: [I,F,T,B]" % (entry.type) )
|
|---|
| 82 | self.header_comments[key] = entry.comment
|
|---|
| 83 |
|
|---|
| 84 | def _setup_columns(self):
|
|---|
| 85 | """
|
|---|
| 86 | """
|
|---|
| 87 | col_type_to_np_type_map = { 'L' : 'b1', 'A' : 'a1', 'B' : 'i1',
|
|---|
| 88 | 'I' : 'i2', 'J' : 'i4', 'K' : 'i8', 'E' : 'f4', 'D' : 'f8'}
|
|---|
| 89 | self.cols = {}
|
|---|
| 90 | for key,col in self.f.GetColumns():
|
|---|
| 91 | self.cols[key] = np.zeros(col.num, col_type_to_np_type_map[col.type])
|
|---|
| 92 | if col.num != 0:
|
|---|
| 93 | self.f.SetPtrAddress(key, self.cols[key])
|
|---|
| 94 |
|
|---|
| 95 | def __iter__(self):
|
|---|
| 96 | return self
|
|---|
| 97 |
|
|---|
| 98 | def next(self, row=None):
|
|---|
| 99 | """
|
|---|
| 100 | """
|
|---|
| 101 | if row is None:
|
|---|
| 102 | if self.f.GetNextRow() == False:
|
|---|
| 103 | raise StopIteration
|
|---|
| 104 | else:
|
|---|
| 105 | row = int(row)
|
|---|
| 106 | if self.f.GetRow(row) == False:
|
|---|
| 107 | raise StopIteration
|
|---|
| 108 | return self
|
|---|
| 109 |
|
|---|
| 110 |
|
|---|
| 111 |
|
|---|
| 112 | class RawData( Fits ):
|
|---|
| 113 | """ Special raw data FITS file access (with DRS4 calibration)
|
|---|
| 114 |
|
|---|
| 115 | During iteration the C++ method DrsCalibration::Apply is being called.
|
|---|
| 116 | """
|
|---|
| 117 | __module__='pyfact'
|
|---|
| 118 | def __init__(self, data_path, calib_path):
|
|---|
| 119 | """ -constructor-
|
|---|
| 120 | *data_path* : fits or fits.gz file of the data including the path
|
|---|
| 121 | *calib_path* : fits or fits.gz file containing DRS calibration data
|
|---|
| 122 | """
|
|---|
| 123 | super(RawData, self).__init__(data_path)
|
|---|
| 124 | self.cols['CalibData'] = np.zeros( self.cols['Data'].shape, np.float32)
|
|---|
| 125 | self.cols['CalibData2D'] = self.cols['CalibData'].reshape( self.header['NPIX'], -1)
|
|---|
| 126 | if not self.cols['CalibData2D'].base is self.cols['CalibData']:
|
|---|
| 127 | print "Error seomthing went wrong!"
|
|---|
| 128 | self.drs_calibration = DrsCalibration()
|
|---|
| 129 | self.drs_calibration.ReadFitsImp( calib_path )
|
|---|
| 130 |
|
|---|
| 131 | self.drs_calibrate = DrsCalibrate()
|
|---|
| 132 | self.list_of_previous_start_cells = []
|
|---|
| 133 |
|
|---|
| 134 | def __iter__(self):
|
|---|
| 135 | """ iterator """
|
|---|
| 136 | return self
|
|---|
| 137 |
|
|---|
| 138 | def next(self, row=None):
|
|---|
| 139 | """
|
|---|
| 140 | """
|
|---|
| 141 | super(RawData, self).next(row)
|
|---|
| 142 |
|
|---|
| 143 | self.drs_calibration.Apply( self.cols['CalibData'],
|
|---|
| 144 | self.cols['Data'],
|
|---|
| 145 | self.cols['StartCellData'],
|
|---|
| 146 | self.header['NROI'])
|
|---|
| 147 |
|
|---|
| 148 | for previous_start_cells in self.list_of_previous_start_cells:
|
|---|
| 149 | self.drs_calibrate.CorrectStep(
|
|---|
| 150 | self.cols['CalibData'],
|
|---|
| 151 | self.header['NPIX'],
|
|---|
| 152 | self.header['NROI'],
|
|---|
| 153 | previous_start_cells,
|
|---|
| 154 | self.cols['StartCellData'],
|
|---|
| 155 | self.header['NROI']+10)
|
|---|
| 156 | self.drs_calibrate.CorrectStep(
|
|---|
| 157 | self.cols['CalibData'],
|
|---|
| 158 | self.header['NPIX'],
|
|---|
| 159 | self.header['NROI'],
|
|---|
| 160 | previous_start_cells,
|
|---|
| 161 | self.cols['StartCellData'],
|
|---|
| 162 | 3)
|
|---|
| 163 | self.list_of_previous_start_cells.append(self.cols['StartCellData'])
|
|---|
| 164 | if len(self.list_of_previous_start_cells) > 5:
|
|---|
| 165 | self.list_of_previous_start_cells.pop(0)
|
|---|
| 166 |
|
|---|
| 167 | for ch in range(self.header['NPIX']):
|
|---|
| 168 | self.drs_calibrate.RemoveSpikes3(self.cols['CalibData2D'][ch], self.header['NROI'])
|
|---|
| 169 |
|
|---|
| 170 | return self
|
|---|
| 171 |
|
|---|
| 172 |
|
|---|
| 173 |
|
|---|
| 174 | if __name__ == '__main__':
|
|---|
| 175 | """ tests """
|
|---|
| 176 | import sys
|
|---|
| 177 | if len(sys.argv) == 2:
|
|---|
| 178 | print "Example for aux-file or uncalibrated raw-file"
|
|---|
| 179 | f = Fits(sys.argv[1])
|
|---|
| 180 | elif len(sys.argv) == 3:
|
|---|
| 181 | print "Example for calibrated raw-file"
|
|---|
| 182 | f = RawData(sys.argv[1], sys.argv[2])
|
|---|
| 183 | else:
|
|---|
| 184 | f = None
|
|---|
| 185 |
|
|---|
| 186 | print "number of events:", f.header['NAXIS2']
|
|---|
| 187 | print "date of observation:", f.header['DATE']
|
|---|
| 188 | print "The files has these cols:", f.cols.keys()
|
|---|
| 189 |
|
|---|
| 190 | for counter,row in enumerate(f):
|
|---|
| 191 | print "Event Id:", row.cols['EventNum']
|
|---|
| 192 | print "shape of column 'StartCellData'", row.cols['StartCellData'].shape
|
|---|
| 193 | print "dtype of column 'Data'", row.cols['StartCellData'].dtype
|
|---|
| 194 | if counter > 3:
|
|---|
| 195 | break
|
|---|
| 196 | # get next row
|
|---|
| 197 | f.next()
|
|---|
| 198 | print "Event Id:", f.cols['EventNum']
|
|---|
| 199 | # get another row
|
|---|
| 200 | f.next(10)
|
|---|
| 201 | print "Event Id:", f.cols['EventNum']
|
|---|
| 202 | # Go back again
|
|---|
| 203 | f.next(3)
|
|---|
| 204 | print "Event Id:", f.cols['EventNum']
|
|---|
| 205 |
|
|---|
| 206 | if len(sys.argv) == 3:
|
|---|
| 207 | import matplotlib.pyplot as plt
|
|---|
| 208 | plt.ion()
|
|---|
| 209 | for i in range(f.header['NPIX']):
|
|---|
| 210 | plt.cla()
|
|---|
| 211 | plt.plot(f.cols['CalibData2D'][i], '.:', label='pixel %d'%i)
|
|---|
| 212 | plt.legend()
|
|---|
| 213 | answer = raw_input('anykey for next pixel; "q" to quit. :')
|
|---|
| 214 | if 'q' in answer.lower():
|
|---|
| 215 | break
|
|---|
| 216 |
|
|---|
