source: fact/tools/pyscripts/new_pyfact/pyfact.py@ 17975

Last change on this file since 17975 was 17975, checked in by dneise, 10 years ago
more pythonic way of getting the number of events
  • Property svn:executable set to *
File size: 9.7 KB
Line 
1#!/usr/bin/python -tt
2#
3# Werner Lustermann, Dominik Neise
4# ETH Zurich, TU Dortmund
5#
6import os
7import sys
8import numpy as np
9import ROOT
10
11########## BUILDING OF THE SHARED OBJECT FILES ###############################
12if __name__ == '__main__' and len(sys.argv) > 1 and 'build' in sys.argv[1]:
13 ROOT.gSystem.AddLinkedLibs("-lz")
14 root_make_string = ROOT.gSystem.GetMakeSharedLib()
15 if not "-std=c++0x" in root_make_string:
16 make_string = root_make_string.replace('$Opt', '$Opt -std=c++0x -D HAVE_ZLIB')
17 ROOT.gSystem.SetMakeSharedLib(make_string)
18 ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/izstream.h+O")
19 ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/fits.h+O")
20 ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/zfits.h+O")
21 ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/factfits.h+O")
22 if not "-std=c++0x" in root_make_string:
23 make_string = root_make_string.replace('$Opt', "$Opt -std=c++0x -D HAVE_ZLIB -D'PACKAGE_NAME=\"PACKAGE_NAME\"' "
24 "-D'PACKAGE_VERSION=\"PACKAGE_VERSION\"' -D'REVISION=\"REVISION\"' ")
25 ROOT.gSystem.SetMakeSharedLib(make_string)
26 ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/DrsCalib.h+O")
27
28 ROOT.gInterpreter.GenerateDictionary("map<string,fits::Entry>","map;string;extern_Mars_mcore/fits.h")
29 ROOT.gInterpreter.GenerateDictionary("pair<string,fits::Entry>","map;string;extern_Mars_mcore/fits.h")
30 ROOT.gInterpreter.GenerateDictionary("map<string,fits::Table::Column>","map;string;extern_Mars_mcore/fits.h")
31 ROOT.gInterpreter.GenerateDictionary("pair<string,fits::Table::Column>","map;string;extern_Mars_mcore/fits.h")
32 ROOT.gInterpreter.GenerateDictionary("vector<DrsCalibrate::Step>","vector;extern_Mars_mcore/DrsCalib.h")
33
34########## USAGE #############################################################
35if __name__ == '__main__' and len(sys.argv) < 3:
36 print """ Usage:
37 ----------------------------------------------------------------------
38 To just build the shared object libs call:
39 python pyfact.py build
40
41 To build (if necessary) and open an example file
42 python -i pyfact.py /path/to/data_file.zfits /path/to/calib_file.drs.fits.gz
43
44 Any of the 3 file 'types': fits.gz zfits fits
45 should be supported.
46
47 To clean all of automatically built files, do something like:
48 rm *.so *.d AutoDict_* extern_Mars_mcore/*.so extern_Mars_mcore/*.d
49 ----------------------------------------------------------------------
50 """
51 sys.exit(1)
52
53######### START OF PYFACT MODULE ############################################
54path = os.path.dirname(os.path.realpath(__file__))
55ROOT.gSystem.Load(path+'/AutoDict_map_string_fits__Entry__cxx.so')
56ROOT.gSystem.Load(path+'/AutoDict_map_string_fits__Table__Column__cxx.so')
57ROOT.gSystem.Load(path+'/AutoDict_pair_string_fits__Entry__cxx.so')
58ROOT.gSystem.Load(path+'/AutoDict_pair_string_fits__Table__Column__cxx.so')
59ROOT.gSystem.Load(path+'/AutoDict_vector_DrsCalibrate__Step__cxx.so')
60ROOT.gSystem.Load(path+'/extern_Mars_mcore/fits_h.so')
61ROOT.gSystem.Load(path+'/extern_Mars_mcore/izstream_h.so')
62ROOT.gSystem.Load(path+'/extern_Mars_mcore/zfits_h.so')
63ROOT.gSystem.Load(path+'/extern_Mars_mcore/factfits_h.so')
64ROOT.gSystem.Load(path+'/extern_Mars_mcore/DrsCalib_h.so')
65del path
66
67class Fits( object ):
68 """ General FITS file access
69
70 Wrapper for factfits class from Mars/mcore/factfits.h
71 (factfits might not be suited to read any FITS file out there, but certainly
72 all FACT FITS files can be read using this class)
73 """
74 __module__ = 'pyfact'
75 def __init__(self, path):
76 """
77 """
78 if not os.path.exists(path):
79 raise IOError(path+' was not found')
80 self.f = ROOT.factfits(path)
81 self._make_header()
82 self._setup_columns()
83
84 def __len__(self):
85 """ return the number of events, in case it's a FACT physics data file
86 or simply the number of rows in case it's a slow data (so called aux) file.
87 """
88 return self.f.GetNumRows()
89
90 def _make_header(self):
91 """
92 """
93 str_to_bool = { 'T':True, 'F':False}
94 type_conversion = { 'I' : int, 'F' : float, 'T' : str, 'B' : str_to_bool.__getitem__}
95
96 self.header = {}
97 self.header_comments = {}
98 for key,entry in self.f.GetKeys():
99 try:
100 self.header[key] = type_conversion[entry.type](entry.value)
101 except KeyError:
102 raise IOError("Error: entry type unknown.\n Is %s, but should be one of: [I,F,T,B]" % (entry.type) )
103 self.header_comments[key] = entry.comment
104
105 def _setup_columns(self):
106 """
107 """
108 col_type_to_np_type_map = { 'L' : 'b1', 'A' : 'a1', 'B' : 'i1',
109 'I' : 'i2', 'J' : 'i4', 'K' : 'i8', 'E' : 'f4', 'D' : 'f8'}
110 self.cols = {}
111 for key,col in self.f.GetColumns():
112 self.cols[key] = np.zeros(col.num, col_type_to_np_type_map[col.type])
113 if col.num != 0:
114 self.f.SetPtrAddress(key, self.cols[key])
115
116 def __iter__(self):
117 return self
118
119 def next(self, row=None):
120 """
121 """
122 if row is None:
123 if self.f.GetNextRow() == False:
124 raise StopIteration
125 else:
126 row = int(row)
127 if self.f.GetRow(row) == False:
128 raise StopIteration
129 return self
130
131class AuxFile( Fits ):
132 """ easy(?) access to FACT aux files
133 """
134 __module__ = 'pyfact'
135 def __init__(self, path, verbose=False):
136 self._verbose = verbose
137 super(AuxFile, self).__init__(path)
138
139 def _setup_columns(self):
140 col_type_to_np_type_map = { 'L' : 'b1', 'A' : 'a1', 'B' : 'i1',
141 'I' : 'i2', 'J' : 'i4', 'K' : 'i8', 'E' : 'f4', 'D' : 'f8'}
142 self._cols = {}
143 self.cols = {}
144
145 N = len(self)
146 for key,col in self.f.GetColumns():
147 self._cols[key] = np.zeros(col.num, col_type_to_np_type_map[col.type])
148 self.cols[key] = np.zeros((N, col.num), col_type_to_np_type_map[col.type])
149 if col.num != 0:
150 self.f.SetPtrAddress(key, self._cols[key])
151
152
153 for i,row in enumerate(self):
154 if self._verbose:
155 try:
156 step = int(self._verbose)
157 except:
158 step = 10
159 if i % step == 0:
160 print "reading line", i, 'of', self.header['NAXIS2']
161
162 for key in self._cols:
163 self.cols[key][i,:] = self._cols[key]
164
165
166
167
168class RawData( Fits ):
169 """ Special raw data FITS file access (with DRS4 calibration)
170
171 During iteration the C++ method DrsCalibration::Apply is being called.
172 """
173 __module__='pyfact'
174 def __init__(self, data_path, calib_path):
175 """ -constructor-
176 *data_path* : fits or fits.gz file of the data including the path
177 *calib_path* : fits or fits.gz file containing DRS calibration data
178 """
179 super(RawData, self).__init__(data_path)
180 self.cols['CalibData'] = np.zeros( self.cols['Data'].shape, np.float32)
181 self.cols['CalibData2D'] = self.cols['CalibData'].reshape( self.header['NPIX'], -1)
182 if not self.cols['CalibData2D'].base is self.cols['CalibData']:
183 print "Error seomthing went wrong!"
184 self.drs_calibration = ROOT.DrsCalibration()
185 self.drs_calibration.ReadFitsImp( calib_path )
186
187 self.drs_calibrate = ROOT.DrsCalibrate()
188 self.list_of_previous_start_cells = []
189
190 def __iter__(self):
191 """ iterator """
192 return self
193
194 def next(self, row=None):
195 """
196 """
197 super(RawData, self).next(row)
198
199 self.drs_calibration.Apply( self.cols['CalibData'],
200 self.cols['Data'],
201 self.cols['StartCellData'],
202 self.header['NROI'])
203
204 for previous_start_cells in self.list_of_previous_start_cells:
205 self.drs_calibrate.CorrectStep(
206 self.cols['CalibData'],
207 self.header['NPIX'],
208 self.header['NROI'],
209 previous_start_cells,
210 self.cols['StartCellData'],
211 self.header['NROI']+10)
212 self.drs_calibrate.CorrectStep(
213 self.cols['CalibData'],
214 self.header['NPIX'],
215 self.header['NROI'],
216 previous_start_cells,
217 self.cols['StartCellData'],
218 3)
219 self.list_of_previous_start_cells.append(self.cols['StartCellData'])
220 if len(self.list_of_previous_start_cells) > 5:
221 self.list_of_previous_start_cells.pop(0)
222
223 for ch in range(self.header['NPIX']):
224 self.drs_calibrate.RemoveSpikes3(self.cols['CalibData2D'][ch], self.header['NROI'])
225
226 return self
227
228
229
230if __name__ == '__main__':
231 """ Example """
232 print "Example for calibrated raw-file"
233 f = RawData(sys.argv[1], sys.argv[2])
234
235 print "number of events:", len(f)
236 print "date of observation:", f.header['DATE']
237 print "The files has these cols:", f.cols.keys()
238
239 for counter,row in enumerate(f):
240 print "Event Id:", row.cols['EventNum']
241 print "shape of column 'StartCellData'", row.cols['StartCellData'].shape
242 print "dtype of column 'Data'", row.cols['StartCellData'].dtype
243 if counter > 3:
244 break
245 # get next row
246 f.next()
247 print "Event Id:", f.cols['EventNum']
248 # get another row
249 f.next(10)
250 print "Event Id:", f.cols['EventNum']
251 # Go back again
252 f.next(3)
253 print "Event Id:", f.cols['EventNum']
Note: See TracBrowser for help on using the repository browser.