Changeset 13124
- Timestamp:
- 03/15/12 22:38:09 (13 years ago)
- File:
-
- 1 edited
-
fact/tools/pyscripts/pyfact/pyfact.py (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
fact/tools/pyscripts/pyfact/pyfact.py
r12947 r13124 4 4 # ETH Zurich 5 5 # 6 from 6 from ctypes import * 7 7 import numpy as np 8 8 from scipy import signal … … 14 14 from ROOT import * 15 15 16 class rawdata( object ): 17 """ 18 raw data access and calibration 16 17 class RawData( object ): 18 """ raw data access and calibration 19 19 20 - open raw data file and drs calibration file 20 21 - performs amplitude calibration … … 23 24 row = number of pixel 24 25 col = length of region of interest 26 25 27 """ 26 # constructor of the classe 27 def __init__( self, dfname, calfname, bslfname='' ): 28 """ 28 29 def __init__(self, data_file_name, calib_file_name, 30 user_action_calib=lambda acal_data, data, blm, tom, gm, scells, nroi: None, 31 baseline_file_name=''): 32 """ initialize object 33 29 34 open data file and calibration data file 30 get basic information about the data in d fname35 get basic information about the data in data_file_name 31 36 allocate buffers for data access 32 37 33 dfname : fits or fits.gz file containing the data including the path 34 calfname : fits or fits.gz file containing DRS calibration data 35 bslfname : npy file containing the baseline values 36 """ 37 self.dfname = dfname 38 self.calfname = calfname 39 self.bslfname = bslfname 38 data_file_name : fits or fits.gz file of the data including the path 39 calib_file_name : fits or fits.gz file containing DRS calibration data 40 baseline_file_name : npy file containing the baseline values 41 42 """ 43 44 self.data_file_name = data_file_name 45 self.calib_file_name = calib_file_name 46 self.baseline_file_name = baseline_file_name 47 48 self.user_action_calib = user_action_calib 40 49 41 50 # baseline correction: True / False 42 if len( bslfname) == 0:51 if len(baseline_file_name) == 0: 43 52 self.correct_baseline = False 44 53 else: … … 47 56 # access data file 48 57 try: 49 d f = fits( self.dfname)58 data_file = fits(self.data_file_name) 50 59 except IOError: 51 print 'problem accessing data file: ', dfname 52 raise # stop ! no data 53 self.df = df 60 print 'problem accessing data file: ', data_file_name 61 raise # stop ! no data 62 #: data file (fits object) 63 self.data_file = data_file 54 64 55 65 # get basic information about the data file 56 self.NROI = df.GetUInt( 'NROI' ) # region of interest (length of DRS pipeline read out) 57 self.NPIX = df.GetUInt( 'NPIX' ) # number of pixels (should be 1440) 58 self.NEvents = df.GetNumRows() # find number of events 66 #: region of interest (number of DRS slices read) 67 self.nroi = data_file.GetUInt('NROI') 68 #: number of pixels (should be 1440) 69 self.npix = data_file.GetUInt('NPIX') 70 #: number of events in the data run 71 self.nevents = data_file.GetNumRows() 72 59 73 # allocate the data memories 60 self.evNum = c_ulong() 61 self.trigType = c_ushort() 62 self.Data = np.zeros( self.NPIX * self.NROI, np.int16 ) 63 self.startCells = np.zeros( self.NPIX, np.int16 ) 74 self.event_id = c_ulong() 75 self.trigger_type = c_ushort() 76 #: 1D array with raw data 77 self.data = np.zeros( self.npix * self.nroi, np.int16 ) 78 #: slice where drs readout started 79 self.start_cells = np.zeros( self.npix, np.int16 ) 80 #: time when the FAD was triggered, in some strange units... 81 self.board_times = np.zeros( 40, np.int32 ) 82 64 83 # set the pointers to the data++ 65 d f.SetPtrAddress( 'EventNum', self.evNum)66 d f.SetPtrAddress( 'TriggerType', self.trigType)67 d f.SetPtrAddress( 'StartCellData', self.startCells ) # DRS readout start cell68 d f.SetPtrAddress( 'Data', self.Data ) # this is what you would expect69 # df.GetNextRow() # access the first event70 71 # accesscalibration file84 data_file.SetPtrAddress('EventNum', self.event_id) 85 data_file.SetPtrAddress('TriggerType', self.trigger_type) 86 data_file.SetPtrAddress('StartCellData', self.start_cells) 87 data_file.SetPtrAddress('Data', self.data) 88 data_file.SetPtrAddress('BoardTime', self.board_times) 89 90 # open the calibration file 72 91 try: 73 cal f = fits( self.calfname)92 calib_file = fits(self.calib_file_name) 74 93 except IOError: 75 print 'problem accessing calibration file: ', cal fname94 print 'problem accessing calibration file: ', calib_file_name 76 95 raise 77 self.calf = calf 78 # 79 BaselineMean = calf.GetN('BaselineMean') 80 GainMean = calf.GetN('GainMean') 81 TriggerOffsetMean = calf.GetN('TriggerOffsetMean') 82 83 self.blm = np.zeros( BaselineMean, np.float32 ) 84 self.gm = np.zeros( GainMean, np.float32 ) 85 self.tom = np.zeros( TriggerOffsetMean, np.float32 ) 86 87 self.Nblm = BaselineMean / self.NPIX 88 self.Ngm = GainMean / self.NPIX 89 self.Ntom = TriggerOffsetMean / self.NPIX 90 91 calf.SetPtrAddress( 'BaselineMean', self.blm ) 92 calf.SetPtrAddress( 'GainMean', self.gm ) 93 calf.SetPtrAddress( 'TriggerOffsetMean', self.tom ) 94 calf.GetRow(0) 95 96 self.v_bsl = np.zeros( self.NPIX ) # array with baseline values (all ZERO) 97 self.smoothData = None 98 self.maxPos = None 99 self.maxAmp = None 100 101 102 def next( self ): 103 """ 104 load the next event from disk and calibrate it 105 """ 106 self.df.GetNextRow() 107 self.calibrate_drsAmplitude() 108 109 110 def calibrate_drsAmplitude( self ): 111 """ 112 perform amplitude calibration for the event 113 """ 114 tomV = 2000./4096. 115 acalData = self.Data * tomV # convert into mV 116 117 # reshape arrays: row = pixel, col = drs_slice 118 acalData = np.reshape( acalData, (self.NPIX, self.NROI) ) 119 blm = np.reshape( self.blm, (self.NPIX, 1024) ) 120 tom = np.reshape( self.tom, (self.NPIX, 1024) ) 121 gm = np.reshape( self.gm, (self.NPIX, 1024) ) 122 123 # print 'acal Data ', acalData.shape 124 # print 'blm shape ', blm.shape 125 # print 'gm shape ', gm.shape 126 127 for pixel in range( self.NPIX ): 96 #: drs calibration file 97 self.calib_file = calib_file 98 99 baseline_mean = calib_file.GetN('BaselineMean') 100 gain_mean = calib_file.GetN('GainMean') 101 trigger_offset_mean = calib_file.GetN('TriggerOffsetMean') 102 103 self.blm = np.zeros(baseline_mean, np.float32) 104 self.gm = np.zeros(gain_mean, np.float32) 105 self.tom = np.zeros(trigger_offset_mean, np.float32) 106 107 self.Nblm = baseline_mean / self.npix 108 self.Ngm = gain_mean / self.npix 109 self.Ntom = trigger_offset_mean / self.npix 110 111 calib_file.SetPtrAddress('BaselineMean', self.blm) 112 calib_file.SetPtrAddress('GainMean', self.gm) 113 calib_file.SetPtrAddress('TriggerOffsetMean', self.tom) 114 calib_file.GetRow(0) 115 116 self.v_bsl = np.zeros(self.npix) # array of baseline values (all ZERO) 117 self.data_saverage_out = None 118 self.pulse_time_of_maximum = None 119 self.pulse_amplitude = None 120 121 def __iter__(self): 122 """ iterator """ 123 return self 124 125 def __add__(self, jump_over): 126 self.data_file.GetRow(jump_over) 127 return self 128 129 def next(self): 130 """ used by __iter__ """ 131 132 if self.data_file.GetNextRow() == False: 133 raise StopIteration 134 else: 135 self.calibrate_drs_amplitude() 136 137 #print 'nevents = ', self.nevents, 'event_id = ', self.event_id.value 138 139 return self.acal_data, self.start_cells, self.trigger_type.value 140 141 def next_event(self): 142 """ load the next event from disk and calibrate it 143 144 """ 145 146 self.data_file.GetNextRow() 147 self.calibrate_drs_amplitude() 148 149 def calibrate_drs_amplitude(self): 150 """ perform the drs amplitude calibration of the event data 151 152 """ 153 154 to_mV = 2000./4096. 155 #: 2D array with amplitude calibrated dat in mV 156 acal_data = self.data * to_mV # convert ADC counts to mV 157 158 # make 2D arrays: row = pixel, col = drs_slice 159 acal_data = np.reshape(acal_data, (self.npix, self.nroi) ) 160 blm = np.reshape(self.blm, (self.npix, self.Nblm) ) 161 tom = np.reshape(self.tom, (self.npix, self.Ntom) ) 162 gm = np.reshape(self.gm, (self.npix, self.Ngm) ) 163 164 for pixel in range( self.npix ): 128 165 # rotate the pixel baseline mean to the Data startCell 129 blm_pixel = np.roll( blm[pixel,:], -self.startCells[pixel] ) 130 acalData[pixel,:] -= blm_pixel[0:self.NROI] 131 acalData[pixel,:] -= tom[pixel, 0:self.NROI] 132 acalData[pixel,:] /= gm[pixel, 0:self.NROI] 166 blm_pixel = np.roll( blm[pixel,:], -self.start_cells[pixel] ) 167 tom_pixel = np.roll( tom[pixel,:], -self.start_cells[pixel] ) 168 gm_pixel = np.roll( gm[pixel,:], -self.start_cells[pixel] ) 169 acal_data[pixel,:] -= blm_pixel[0:self.nroi] 170 acal_data[pixel,:] -= tom_pixel[0:self.nroi] 171 acal_data[pixel,:] /= gm_pixel[0:self.nroi] 133 172 134 self.acalData = acalData * 1907.35 135 136 # print 'acalData ', self.acalData[0:2,0:20] 137 138 def filterSlidingAverage( self , windowSize = 4): 139 """ sliding average filter 140 using: 141 self.acalData 142 filling array: 143 self.smoothData 144 """ 145 #scipy.signal.lfilter(b, a, x, axis=-1, zi=None) 146 smoothData = self.acalData.copy() 147 b = np.ones( windowSize ) 148 a = np.zeros( windowSize ) 149 a[0] = len(b) 150 smoothData[:,:] = signal.lfilter(b, a, smoothData[:,:]) 151 152 self.smoothData = smoothData 153 154 def filterCFD( self, length=10, ratio=0.75): 155 """ constant fraction filter 156 using: 157 self.smoothData 158 filling array: 159 self.cfdData 160 """ 161 if self.smoothData == None: 162 print 'error pyfact.filterCFD was called without prior call to filterSlidingAverage' 163 print ' variable self.smoothData is needed ' 164 pass 165 cfdData = self.smoothData.copy() 166 b = np.zeros( length ) 167 a = np.zeros( length ) 168 b[0] = -1. * ratio 169 b[length-1] = 1. 170 a[0] = 1. 171 cfdData[:,:] = signal.lfilter(b, a, cfdData[:,:]) 172 173 self.cfdData = cfdData 174 175 def findPeak (self, min=30, max=250): 176 """ find maximum in search window 177 using: 178 self.smoothData 179 filling arrays: 180 self.maxPos 181 self.maxAmp 182 """ 183 if self.smoothData == None: 184 print 'error pyfact.findPeakMax was called without prior call to filterSlidingAverage' 185 print ' variable self.smoothData is needed ' 186 pass 187 maxPos = np.argmax( self.smoothData[:,min:max] , 1) 188 maxAmp = np.max( self.smoothData[:,min:max] , 1) 189 self.maxPos = maxPos 190 self.maxAmp = maxAmp 191 192 def sumAroundPeak (self, left=13, right=23): 193 """ integrate signal in gate around Peak 194 using: 195 self.maxPos 196 self.acalData 197 filling array: 198 self.integral 199 """ 200 if self.maxPos == None: 201 print 'error pyfact.sumAroundPeak was called without prior call of findPeak' 202 print ' variable self.maxPos is needed' 203 pass 204 205 sums = np.empty( self.NPIX ) 206 for pixel in range( self.NPIX ): 207 min = self.maxPos[pixel]-left 208 max = self.maxPos[pixel]+right 209 sums[pixel] = self.acalData[pixel,min:max].sum() 210 211 self.integral = sums 212 213 def ReadBaseline( self, file, bsl_hist = 'bsl_sum/hplt_mean' ): 214 """ 173 self.acal_data = acal_data * 1907.35 174 175 #print 'blm _pyfact', blm[0,0:20] 176 #t = np.roll( blm[0,:], -self.start_cells[0] ) 177 #print 'blm _pyfact', t[0:20] 178 #print 'start_pyfact: ', self.start_cells[0] 179 #print 'acal _pyfact: ', self.acal_data[0,0:10] 180 #t = np.roll( gm[0,:], -self.start_cells[0] ) 181 #print 'gm _pyfact: ', t[0:10] 182 self.user_action_calib( self.acal_data, 183 np.reshape(self.data, (self.npix, self.nroi) ), blm, tom, gm, self.start_cells, self.nroi) 184 185 186 def baseline_read_values(self, file, bsl_hist='bsl_sum/hplt_mean'): 187 """ 188 215 189 open ROOT file with baseline histogram and read baseline values 216 190 file name of the root file 217 191 bsl_hist path to the histogram containing the basline values 218 """ 192 193 """ 194 219 195 try: 220 f = TFile( file)196 f = TFile(file) 221 197 except: 222 198 print 'Baseline data file could not be read: ', file 223 199 return 224 200 225 h = f.Get( bsl_hist)226 227 for i in range( self.NPIX):228 self.v_bsl[i] = h.GetBinContent( i+1)201 h = f.Get(bsl_hist) 202 203 for i in range(self.npix): 204 self.v_bsl[i] = h.GetBinContent(i+1) 229 205 230 206 f.Close() 231 232 233 def CorrectBaseline( self ):234 """ 235 apply baseline correction236 """237 for pixel in range( self.NPIX):238 self.acal Data[pixel,:] -= self.v_bsl[pixel]239 240 241 def info( self ):242 """243 print information244 """245 print 'data file: ', d fname246 print 'calib file: ', cal fname207 208 def baseline_correct(self): 209 """ subtract baseline from the data 210 211 """ 212 213 for pixel in range(self.npix): 214 self.acal_data[pixel,:] -= self.v_bsl[pixel] 215 216 def info(self): 217 """ print run information 218 219 """ 220 221 print 'data file: ', data_file_name 222 print 'calib file: ', calib_file_name 247 223 print 'calibration file' 248 print 'N BaselineMean: ', self.Nblm249 print 'N GainMean: ', self.Ngm224 print 'N baseline_mean: ', self.Nblm 225 print 'N gain mean: ', self.Ngm 250 226 print 'N TriggeroffsetMean: ', self.Ntom 251 227 252 # ----------------------------------------------------------------------------- ---228 # ----------------------------------------------------------------------------- 253 229 class fnames( object ): 254 230 """ organize file names of a FACT data run … … 256 232 """ 257 233 258 def __init__( 234 def __init__(self, specifier = ['012', '023', '2011', '11', '24'], 259 235 rpath = '/scratch_nfs/res/bsl/', 260 236 zipped = True): … … 265 241 rpath : directory path for the results; YYYYMMDD will be appended to rpath 266 242 zipped : use zipped (True) or unzipped (Data) 267 """ 243 244 """ 245 268 246 self.specifier = specifier 269 247 self.rpath = rpath … … 271 249 272 250 self.make( self.specifier, self.rpath, self.zipped ) 273 # end of def __init__ 251 274 252 275 253 def make( self, specifier, rpath, zipped ): … … 312 290 self.results = self.names['results'] 313 291 314 # end of make315 316 292 def info( self ): 317 293 """ print complete filenames … … 323 299 print 'drs-cal: ', self.names['drscal'] 324 300 print 'results: ', self.names['results'] 325 # end of def info326 301 327 302 # end of class definition: fnames( object ) 328 303 329 330 331 class histogramList( object ): 332 333 def __init__( self, name ): 334 """ set the name and create empty lists """ 335 self.name = name # name of the list 336 self.list = [] # list of the histograms 337 self.dict = {} # dictionary of histograms 338 self.hList = TObjArray() # list a la ROOT of the histograms 339 340 def add( self, tag, h ): 341 self.list.append( h ) 342 self.dict[tag] = h 343 self.hList.Add( h ) 344 345 346 class pixelHisto1d ( object ): 347 348 def __init__( self, name, title, Nbin, first, last, xtitle, ytitle, NPIX ): 349 """ 350 book one dimensional histograms for each pixel 351 """ 352 self.name = name 353 354 self.list = [ x for x in range( NPIX ) ] 355 self.hList = TObjArray() 356 357 for pixel in range( NPIX ): 358 359 hname = name + ' ' + str( pixel ) 360 htitle = title + ' ' + str( pixel ) 361 self.list[pixel] = TH1F( hname, htitle, Nbin, first, last ) 362 363 self.list[pixel].GetXaxis().SetTitle( xtitle ) 364 self.list[pixel].GetYaxis().SetTitle( ytitle ) 365 self.hList.Add( self.list[pixel] ) 366 367 # simple test method 304 def _test_iter(): 305 """ test for function __iter__ """ 306 307 # data_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_111.fits.gz' 308 # calib_file_name = data_file_name = 309 data_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_028.fits.gz' 310 calib_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_022.drs.fits.gz' 311 run = RawData( data_file_name, calib_file_name ) 312 313 for data, scell, tt in run: 314 print 'data[0,0] = ', data[0,0], "start_cell[0] =", scell[0], "trigger type = ", tt 315 316 368 317 if __name__ == '__main__': 369 """ 370 create an instance 371 """ 372 dfname = '/data03/fact-construction/raw/2011/11/24/20111124_121.fits' 373 calfname = '/data03/fact-construction/raw/2011/11/24/20111124_111.drs.fits' 374 rd = rawdata( dfname, calfname ) 375 rd.info() 376 rd.next() 377 378 # for i in range(10): 379 # df.GetNextRow() 380 381 # print 'evNum: ', evNum.value 382 # print 'startCells[0:9]: ', startCells[0:9] 383 # print 'evData[0:9]: ', evData[0:9] 318 """ tests """ 319 320 _test_iter()
Note:
See TracChangeset
for help on using the changeset viewer.
