| 1 | #!/usr/bin/python
|
|---|
| 2 | #
|
|---|
| 3 | # Werner Lustermann
|
|---|
| 4 | # ETH Zurich
|
|---|
| 5 | #
|
|---|
| 6 | from ctypes import *
|
|---|
| 7 | import numpy as np
|
|---|
| 8 |
|
|---|
| 9 | # get the ROOT stuff + my shared libs
|
|---|
| 10 | from ROOT import gSystem
|
|---|
| 11 | # fitslib.so is made from fits.h and is used to access the data
|
|---|
| 12 | gSystem.Load('~/py/fitslib.so')
|
|---|
| 13 | from ROOT import *
|
|---|
| 14 |
|
|---|
| 15 | class rawdata( object ):
|
|---|
| 16 | """
|
|---|
| 17 | raw data access and calibration
|
|---|
| 18 | - open raw data file and drs calibration file
|
|---|
| 19 | - performs amplitude calibration
|
|---|
| 20 | - performs baseline substraction if wanted
|
|---|
| 21 | - provides all data in an array:
|
|---|
| 22 | row = number of pixel
|
|---|
| 23 | col = length of region of interest
|
|---|
| 24 | """
|
|---|
| 25 | # constructor of the classe
|
|---|
| 26 | def __init__( self, dfname, calfname, bslfname='' ):
|
|---|
| 27 | """
|
|---|
| 28 | open data file and calibration data file
|
|---|
| 29 | get basic information about the data in dfname
|
|---|
| 30 | allocate buffers for data access
|
|---|
| 31 |
|
|---|
| 32 | dfname : fits or fits.gz file containing the data including the path
|
|---|
| 33 | calfname : fits or fits.gz file containing DRS calibration data
|
|---|
| 34 | bslfname : npy file containing the baseline values
|
|---|
| 35 | """
|
|---|
| 36 | self.dfname = dfname
|
|---|
| 37 | self.calfname = calfname
|
|---|
| 38 | self.bslfname = bslfname
|
|---|
| 39 |
|
|---|
| 40 | # baseline correction: True / False
|
|---|
| 41 | if len( bslfname ) == 0:
|
|---|
| 42 | self.correct_baseline = False
|
|---|
| 43 | else:
|
|---|
| 44 | self.correct_baseline = True
|
|---|
| 45 |
|
|---|
| 46 | # access data file
|
|---|
| 47 | try:
|
|---|
| 48 | df = fits( self.dfname )
|
|---|
| 49 | except IOError:
|
|---|
| 50 | print 'problem accessing data file: ', dfname
|
|---|
| 51 | raise # stop ! no data
|
|---|
| 52 | self.df = df
|
|---|
| 53 |
|
|---|
| 54 | # get basic information about the data file
|
|---|
| 55 | self.NROI = df.GetUInt( 'NROI' ) # region of interest (length of DRS pipeline read out)
|
|---|
| 56 | self.NPIX = df.GetUInt( 'NPIX' ) # number of pixels (should be 1440)
|
|---|
| 57 | self.NEvents = df.GetNumRows() # find number of events
|
|---|
| 58 | # allocate the data memories
|
|---|
| 59 | self.evNum = c_ulong()
|
|---|
| 60 | self.Data = np.zeros( self.NPIX * self.NROI, np.int16 )
|
|---|
| 61 | self.startCells = np.zeros( self.NPIX, np.int16 )
|
|---|
| 62 | # set the pointers to the data++
|
|---|
| 63 | df.SetPtrAddress( 'EventNum', self.evNum )
|
|---|
| 64 | df.SetPtrAddress( 'StartCellData', self.startCells ) # DRS readout start cell
|
|---|
| 65 | df.SetPtrAddress( 'Data', self.Data ) # this is what you would expect
|
|---|
| 66 | # df.GetNextRow() # access the first event
|
|---|
| 67 |
|
|---|
| 68 | # access calibration file
|
|---|
| 69 | try:
|
|---|
| 70 | calf = fits( self.calfname )
|
|---|
| 71 | except IOError:
|
|---|
| 72 | print 'problem accessing calibration file: ', calfname
|
|---|
| 73 | raise
|
|---|
| 74 | self.calf = calf
|
|---|
| 75 | #
|
|---|
| 76 | BaselineMean = calf.GetN('BaselineMean')
|
|---|
| 77 | GainMean = calf.GetN('GainMean')
|
|---|
| 78 | TriggerOffsetMean = calf.GetN('TriggerOffsetMean')
|
|---|
| 79 |
|
|---|
| 80 | self.blm = np.zeros( BaselineMean, np.float32 )
|
|---|
| 81 | self.gm = np.zeros( GainMean, np.float32 )
|
|---|
| 82 | self.tom = np.zeros( TriggerOffsetMean, np.float32 )
|
|---|
| 83 |
|
|---|
| 84 | self.Nblm = BaselineMean / self.NPIX
|
|---|
| 85 | self.Ngm = GainMean / self.NPIX
|
|---|
| 86 | self.Ntom = TriggerOffsetMean / self.NPIX
|
|---|
| 87 |
|
|---|
| 88 | calf.SetPtrAddress( 'BaselineMean', self.blm )
|
|---|
| 89 | calf.SetPtrAddress( 'GainMean', self.gm )
|
|---|
| 90 | calf.SetPtrAddress( 'TriggerOffsetMean', self.tom )
|
|---|
| 91 | calf.GetRow(0)
|
|---|
| 92 |
|
|---|
| 93 | self.v_bsl = np.zeros( self.NPIX ) # array with baseline values (all ZERO)
|
|---|
| 94 |
|
|---|
| 95 |
|
|---|
| 96 | def next( self ):
|
|---|
| 97 | """
|
|---|
| 98 | load the next event from disk and calibrate it
|
|---|
| 99 | """
|
|---|
| 100 | self.df.GetNextRow()
|
|---|
| 101 | self.calibrate_drsAmplitude()
|
|---|
| 102 |
|
|---|
| 103 |
|
|---|
| 104 | def calibrate_drsAmplitude( self ):
|
|---|
| 105 | """
|
|---|
| 106 | perform amplitude calibration for the event
|
|---|
| 107 | """
|
|---|
| 108 | tomV = 2000./4096.
|
|---|
| 109 | acalData = self.Data * tomV # convert into mV
|
|---|
| 110 |
|
|---|
| 111 | # reshape arrays: row = pixel, col = drs_slice
|
|---|
| 112 | acalData = np.reshape( acalData, (self.NPIX, self.NROI) )
|
|---|
| 113 | blm = np.reshape( self.blm, (self.NPIX, 1024) )
|
|---|
| 114 | tom = np.reshape( self.tom, (self.NPIX, 1024) )
|
|---|
| 115 | gm = np.reshape( self.gm, (self.NPIX, 1024) )
|
|---|
| 116 |
|
|---|
| 117 | # print 'acal Data ', acalData.shape
|
|---|
| 118 | # print 'blm shape ', blm.shape
|
|---|
| 119 | # print 'gm shape ', gm.shape
|
|---|
| 120 |
|
|---|
| 121 | for pixel in range( self.NPIX ):
|
|---|
| 122 | # rotate the pixel baseline mean to the Data startCell
|
|---|
| 123 | blm_pixel = np.roll( blm[pixel,:], -self.startCells[pixel] )
|
|---|
| 124 | acalData[pixel,:] -= blm_pixel[0:self.NROI]
|
|---|
| 125 | acalData[pixel,:] -= tom[pixel, 0:self.NROI]
|
|---|
| 126 | acalData[pixel,:] /= gm[pixel, 0:self.NROI]
|
|---|
| 127 |
|
|---|
| 128 | self.acalData = acalData * 1907.35
|
|---|
| 129 |
|
|---|
| 130 | # print 'acalData ', self.acalData[0:2,0:20]
|
|---|
| 131 |
|
|---|
| 132 | def ReadBaseline( self, file, bsl_hist = 'bsl_sum/hplt_mean' ):
|
|---|
| 133 | """
|
|---|
| 134 | open ROOT file with baseline histogram and read baseline values
|
|---|
| 135 | file name of the root file
|
|---|
| 136 | bsl_hist path to the histogram containing the basline values
|
|---|
| 137 | """
|
|---|
| 138 | try:
|
|---|
| 139 | f = TFile( file )
|
|---|
| 140 | except:
|
|---|
| 141 | print 'Baseline data file could not be read: ', file
|
|---|
| 142 | return
|
|---|
| 143 |
|
|---|
| 144 | h = f.Get( bsl_hist )
|
|---|
| 145 |
|
|---|
| 146 | for i in range( self.NPIX ):
|
|---|
| 147 | self.v_bsl[i] = h.GetBinContent( i+1 )
|
|---|
| 148 |
|
|---|
| 149 | f.Close()
|
|---|
| 150 |
|
|---|
| 151 |
|
|---|
| 152 | def CorrectBaseline( self ):
|
|---|
| 153 | """
|
|---|
| 154 | apply baseline correction
|
|---|
| 155 | """
|
|---|
| 156 | for pixel in range( self.NPIX ):
|
|---|
| 157 | self.acalData[pixel,:] -= self.v_bsl[pixel]
|
|---|
| 158 |
|
|---|
| 159 |
|
|---|
| 160 | def info( self ):
|
|---|
| 161 | """
|
|---|
| 162 | print information
|
|---|
| 163 | """
|
|---|
| 164 | print 'data file: ', dfname
|
|---|
| 165 | print 'calib file: ', calfname
|
|---|
| 166 | print 'calibration file'
|
|---|
| 167 | print 'N BaselineMean: ', self.Nblm
|
|---|
| 168 | print 'N GainMean: ', self.Ngm
|
|---|
| 169 | print 'N TriggeroffsetMean: ', self.Ntom
|
|---|
| 170 |
|
|---|
| 171 | # --------------------------------------------------------------------------------
|
|---|
| 172 | class fnames( object ):
|
|---|
| 173 | """ organize file names of a FACT data run
|
|---|
| 174 |
|
|---|
| 175 | """
|
|---|
| 176 |
|
|---|
| 177 | def __init__( self, specifier = ['012', '023', '2011', '11', '24'],
|
|---|
| 178 | rpath = '/scratch_nfs/bsl/',
|
|---|
| 179 | zipped = True):
|
|---|
| 180 | """
|
|---|
| 181 | specifier : list of strings defined as:
|
|---|
| 182 | [ 'DRS calibration file', 'Data file', 'YYYY', 'MM', 'DD']
|
|---|
| 183 |
|
|---|
| 184 | rpath : directory path for the results; YYYYMMDD will be appended to rpath
|
|---|
| 185 | zipped : use zipped (True) or unzipped (Data)
|
|---|
| 186 | """
|
|---|
| 187 | self.specifier = specifier
|
|---|
| 188 | self.rpath = rpath
|
|---|
| 189 | self.zipped = zipped
|
|---|
| 190 |
|
|---|
| 191 | self.make( self.specifier, self.rpath, self.zipped )
|
|---|
| 192 | # end of def __init__
|
|---|
| 193 |
|
|---|
| 194 | def make( self, specifier, rpath, zipped ):
|
|---|
| 195 | """ create (make) the filenames
|
|---|
| 196 |
|
|---|
| 197 | names : dictionary of filenames, tags { 'data', 'drscal', 'results' }
|
|---|
| 198 | data : name of the data file
|
|---|
| 199 | drscal : name of the drs calibration file
|
|---|
| 200 | results : radikal of file name(s) for results (to be completed by suffixes)
|
|---|
| 201 | """
|
|---|
| 202 |
|
|---|
| 203 | self.specifier = specifier
|
|---|
| 204 |
|
|---|
| 205 | if zipped:
|
|---|
| 206 | dpath = '/data00/fact-construction/raw/'
|
|---|
| 207 | ext = '.fits.gz'
|
|---|
| 208 | else:
|
|---|
| 209 | dpath = '/data03/fact-construction/raw/'
|
|---|
| 210 | ext = '.fits'
|
|---|
| 211 |
|
|---|
| 212 | year = specifier[2]
|
|---|
| 213 | month = specifier[3]
|
|---|
| 214 | day = specifier[4]
|
|---|
| 215 |
|
|---|
| 216 | yyyymmdd = year + month + day
|
|---|
| 217 | dfile = specifier[1]
|
|---|
| 218 | cfile = specifier[0]
|
|---|
| 219 |
|
|---|
| 220 | rpath = rpath + yyyymmdd + '/'
|
|---|
| 221 | self.rpath = rpath
|
|---|
| 222 | self.names = {}
|
|---|
| 223 |
|
|---|
| 224 | tmp = dpath + year + '/' + month + '/' + day + '/' + yyyymmdd + '_'
|
|---|
| 225 | self.names['data'] = tmp + dfile + ext
|
|---|
| 226 | self.names['drscal'] = tmp + '_' + cfile + '.drs' + ext
|
|---|
| 227 | self.names['results'] = rpath + yyyymmdd + '_' + dfile + '_' + cfile
|
|---|
| 228 |
|
|---|
| 229 | self.data = self.names['data']
|
|---|
| 230 | self.drscal = self.names['drscal']
|
|---|
| 231 | self.results = self.names['results']
|
|---|
| 232 |
|
|---|
| 233 | # end of make
|
|---|
| 234 |
|
|---|
| 235 | def info( self ):
|
|---|
| 236 | """ print complete filenames
|
|---|
| 237 |
|
|---|
| 238 | """
|
|---|
| 239 |
|
|---|
| 240 | print 'file names:'
|
|---|
| 241 | print 'data: ', self.names['data']
|
|---|
| 242 | print 'drs-cal: ', self.names['drscal']
|
|---|
| 243 | print 'results: ', self.names['results']
|
|---|
| 244 | # end of def info
|
|---|
| 245 |
|
|---|
| 246 | # end of class definition: fnames( object )
|
|---|
| 247 |
|
|---|
| 248 |
|
|---|
| 249 |
|
|---|
| 250 | class histogramList( object ):
|
|---|
| 251 |
|
|---|
| 252 | def __init__( self, name ):
|
|---|
| 253 | """ set the name and create empty lists """
|
|---|
| 254 | self.name = name # name of the list
|
|---|
| 255 | self.list = [] # list of the histograms
|
|---|
| 256 | self.dict = {} # dictionary of histograms
|
|---|
| 257 | self.hList = TObjArray() # list a la ROOT of the histograms
|
|---|
| 258 |
|
|---|
| 259 | def add( self, tag, h ):
|
|---|
| 260 | self.list.append( h )
|
|---|
| 261 | self.dict[tag] = h
|
|---|
| 262 | self.hList.Add( h )
|
|---|
| 263 |
|
|---|
| 264 |
|
|---|
| 265 | class pixelHisto1d ( object ):
|
|---|
| 266 |
|
|---|
| 267 | def __init__( self, name, title, Nbin, first, last, xtitle, ytitle, NPIX ):
|
|---|
| 268 | """
|
|---|
| 269 | book one dimensional histograms for each pixel
|
|---|
| 270 | """
|
|---|
| 271 | self.name = name
|
|---|
| 272 |
|
|---|
| 273 | self.list = [ x for x in range( NPIX ) ]
|
|---|
| 274 | self.hList = TObjArray()
|
|---|
| 275 |
|
|---|
| 276 | for pixel in range( NPIX ):
|
|---|
| 277 |
|
|---|
| 278 | hname = name + ' ' + str( pixel )
|
|---|
| 279 | htitle = title + ' ' + str( pixel )
|
|---|
| 280 | self.list[pixel] = TH1F( hname, htitle, Nbin, first, last )
|
|---|
| 281 |
|
|---|
| 282 | self.list[pixel].GetXaxis().SetTitle( xtitle )
|
|---|
| 283 | self.list[pixel].GetYaxis().SetTitle( ytitle )
|
|---|
| 284 | self.hList.Add( self.list[pixel] )
|
|---|
| 285 |
|
|---|
| 286 | # simple test method
|
|---|
| 287 | if __name__ == '__main__':
|
|---|
| 288 | """
|
|---|
| 289 | create an instance
|
|---|
| 290 | """
|
|---|
| 291 | dfname = '/data03/fact-construction/raw/2011/11/24/20111124_121.fits'
|
|---|
| 292 | calfname = '/data03/fact-construction/raw/2011/11/24/20111124_111.drs.fits'
|
|---|
| 293 | rd = rawdata( dfname, calfname )
|
|---|
| 294 | rd.info()
|
|---|
| 295 | rd.next()
|
|---|
| 296 |
|
|---|
| 297 | # for i in range(10):
|
|---|
| 298 | # df.GetNextRow()
|
|---|
| 299 |
|
|---|
| 300 | # print 'evNum: ', evNum.value
|
|---|
| 301 | # print 'startCells[0:9]: ', startCells[0:9]
|
|---|
| 302 | # print 'evData[0:9]: ', evData[0:9]
|
|---|
