Index: /fact/tools/pyscripts/pyfact.py
===================================================================
--- /fact/tools/pyscripts/pyfact.py	(revision 13245)
+++ /fact/tools/pyscripts/pyfact.py	(revision 13245)
@@ -0,0 +1,400 @@
+#!/usr/bin/python
+#
+# Werner Lustermann
+# ETH Zurich
+#
+from ctypes import *
+import numpy as np
+from scipy import signal
+
+# get the ROOT stuff + my shared libs
+from ROOT import gSystem
+# fitslib.so is made from fits.h and is used to access the data
+gSystem.Load('~/py/fitslib.so')
+from ROOT import *
+
+
+class RawData( object ):
+    """ raw data access and calibration
+        - open raw data file and drs calibration file
+        - performs amplitude calibration
+        - performs baseline substraction if wanted
+        - provides all data in an array:
+          row = number of pixel
+          col = length of region of interest
+    """
+
+    def __init__(self, data_file_name, calib_file_name, 
+                user_action_calib=lambda acal_data, data, blm, tom, gm, scells, nroi: None,
+                baseline_file_name='',
+                return_dict = False):
+        """ initialize object
+            open data file and calibration data file
+            get basic information about the data in data_file_name
+            allocate buffers for data access
+
+            data_file_name   : fits or fits.gz file of the data including the path
+            calib_file_name : fits or fits.gz file containing DRS calibration data
+            baseline_file_name : npy file containing the baseline values
+        """
+        self.return_dict = return_dict
+        if return_dict == False:
+            print 'Warning: Rawdata.__init__() default value of return_dict is False.'
+            print 'THIS DEFAULT VALUE WILL CHANGE TO **TRUE** SOON.'
+            print
+            print 'After a while the option will be deleted, please change all'
+            print 'your scripts such they use the dictionary output.'
+        self.info = {}      # dictionary containing the entire information 
+        # shortcut
+        info = self.info
+        
+        self.data_file_name = data_file_name
+        info['data_file_name'] = self.data_file_name
+        self.calib_file_name = calib_file_name
+        info['calib_file_name'] = self.calib_file_name
+        self.baseline_file_name = baseline_file_name
+        info['baseline_file_name'] = self.baseline_file_name
+
+        self.user_action_calib = user_action_calib
+        info['user_action_calib'] = self.user_action_calib
+        
+        # baseline correction: True / False
+        if len(baseline_file_name) == 0:
+            self.correct_baseline = False
+        else:
+            self.correct_baseline = True
+        info['correct_baseline'] = self.correct_baseline
+        
+        # access data file
+        try:
+            data_file = fits(self.data_file_name)
+        except IOError:
+            print 'problem accessing data file: ', data_file_name
+            raise  # stop ! no data
+        #: data file (fits object)
+        self.data_file = data_file
+        
+        # get basic information about the data file
+        #: region of interest (number of DRS slices read)
+        self.nroi    = data_file.GetUInt('NROI')
+        info['nroi'] = self.nroi
+        #: number of pixels (should be 1440)
+        self.npix    = data_file.GetUInt('NPIX')
+        info['npix'] = self.npix
+        #: number of events in the data run
+        self.nevents = data_file.GetNumRows()
+        info['nevents'] = self.nevents
+        
+        # allocate the data memories
+        self.event_id = c_ulong()
+        info.['event_id'] = self.event_id.value
+        self.trigger_type = c_ushort()
+        info['trigger_type'] = self.trigger_type.value
+        #: 1D array with raw data
+        self.data  = np.zeros( self.npix * self.nroi, np.int16 )
+        info['raw_data'] = self.data
+        #: slice where drs readout started
+        self.start_cells = np.zeros( self.npix, np.int16 )
+        info['start_cells'] = self.start_cells
+        #: time when the FAD was triggered, in some strange units...
+        self.board_times = np.zeros( 40, np.int32 )
+        info['board_times'] = self.board_times
+
+        # set the pointers to the data++
+        data_file.SetPtrAddress('EventNum', self.event_id)
+        data_file.SetPtrAddress('TriggerType', self.trigger_type)
+        data_file.SetPtrAddress('StartCellData', self.start_cells) 
+        data_file.SetPtrAddress('Data', self.data) 
+        data_file.SetPtrAddress('BoardTime', self.board_times) 
+                
+        # open the calibration file
+        try:
+            calib_file = fits(self.calib_file_name)
+        except IOError:
+            print 'problem accessing calibration file: ', calib_file_name
+            raise
+        #: drs calibration file
+        self.calib_file = calib_file
+        
+        baseline_mean = calib_file.GetN('BaselineMean')
+        gain_mean = calib_file.GetN('GainMean')
+        trigger_offset_mean = calib_file.GetN('TriggerOffsetMean')
+
+        self.blm = np.zeros(baseline_mean, np.float32)
+        self.gm  = np.zeros(gain_mean, np.float32)
+        self.tom = np.zeros(trigger_offset_mean, np.float32)
+        info['blm'] = self.blm
+        info['gm'] = self.gm
+        info['tom'] = self.tom
+
+        self.Nblm = baseline_mean / self.npix
+        self.Ngm  = gain_mean / self.npix
+        self.Ntom  = trigger_offset_mean / self.npix
+        info['Nblm'] = self.Nblm
+        info['Ngm'] = self.Ngm
+        info['Ntom'] = self.Ntom
+
+        calib_file.SetPtrAddress('BaselineMean', self.blm)
+        calib_file.SetPtrAddress('GainMean', self.gm)
+        calib_file.SetPtrAddress('TriggerOffsetMean', self.tom)
+        calib_file.GetRow(0)
+
+        self.v_bsl = np.zeros(self.npix)  # array of baseline values (all ZERO)
+        info['v_bsl'] = self.v_bsl
+
+    def __iter__(self):
+        """ iterator """
+        return self
+        
+    def __add__(self, jump_over):
+        self.data_file.GetRow(jump_over)
+        return self
+        
+    def next(self):
+        """ used by __iter__ """
+
+        if self.data_file.GetNextRow() == False:
+            raise StopIteration
+        else:
+            self.calibrate_drs_amplitude()
+
+        #print 'nevents = ', self.nevents, 'event_id = ', self.event_id.value
+        if self.return_dict:
+            return self.info
+        else:
+            return self.acal_data, self.start_cells, self.trigger_type.value
+
+    def next_event(self):
+        """ load the next event from disk and calibrate it
+        
+        """
+
+        self.data_file.GetNextRow()
+        self.calibrate_drs_amplitude()
+
+    def calibrate_drs_amplitude(self):
+        """ perform the drs amplitude calibration of the event data
+        
+        """
+        # shortcut
+        info = self.info
+
+        to_mV = 2000./4096.
+        #: 2D array with amplitude calibrated dat in mV
+        acal_data = self.data * to_mV  # convert ADC counts to mV
+
+        # make 2D arrays: row = pixel, col = drs_slice
+        acal_data = np.reshape(acal_data, (self.npix, self.nroi) )
+        blm = np.reshape(self.blm, (self.npix, self.Nblm) )
+        tom = np.reshape(self.tom, (self.npix, self.Ntom) )
+        gm  = np.reshape(self.gm,  (self.npix, self.Ngm) )
+        
+        for pixel in range( self.npix ):
+            # rotate the pixel baseline mean to the Data startCell
+            blm_pixel = np.roll( blm[pixel,:], -self.start_cells[pixel] )
+            # rotate the pixel gain mean to the Data startCell
+            gm_pixel = np.roll( gm[pixel,:], -self.start_cells[pixel] )
+            # the 'trigger offset mean' does not need to be rolled
+            # on the contrary, it seems there is an offset in the DRS data,
+            # which is related to its distance to the startCell, not to its 
+            # distance to the beginning of the physical pipeline in the DRS chip
+            tom_pixel = tom[pixel,:]
+            
+            acal_data[pixel,:] -= blm_pixel[0:self.nroi]
+            acal_data[pixel,:] -= tom_pixel[0:self.nroi]
+            acal_data[pixel,:] /= gm_pixel[0:self.nroi]
+            
+        self.acal_data = acal_data * 1907.35
+        info['data'] = self.acal_data
+        #print 'blm _pyfact', blm[0,0:20] 
+        #t = np.roll( blm[0,:], -self.start_cells[0] )
+        #print 'blm _pyfact', t[0:20]
+        #print 'start_pyfact: ', self.start_cells[0]
+        #print 'acal _pyfact: ', self.acal_data[0,0:10] 
+        #t = np.roll( gm[0,:], -self.start_cells[0] )
+        #print 'gm _pyfact: ', t[0:10] 
+        self.user_action_calib( self.acal_data, 
+            np.reshape(self.data, (self.npix, self.nroi) ), blm, tom, gm, self.start_cells, self.nroi)
+
+        
+    def baseline_read_values(self, file, bsl_hist='bsl_sum/hplt_mean'):
+        """
+        
+        open ROOT file with baseline histogram and read baseline values
+        file       name of the root file
+        bsl_hist   path to the histogram containing the basline values
+
+        """
+
+        try:
+            f = TFile(file)
+        except:
+            print 'Baseline data file could not be read: ', file
+            return
+        
+        h = f.Get(bsl_hist)
+
+        for i in range(self.npix):
+            self.v_bsl[i] = h.GetBinContent(i+1)
+
+        f.Close()
+        
+    def baseline_correct(self):
+        """ subtract baseline from the data
+
+        """
+        
+        for pixel in range(self.npix):
+            self.acal_data[pixel,:] -= self.v_bsl[pixel]
+                    
+    def info(self):
+        """ print run information
+        
+        """
+        
+        print 'data file:  ', data_file_name
+        print 'calib file: ', calib_file_name
+        print 'calibration file'
+        print 'N baseline_mean: ', self.Nblm
+        print 'N gain mean: ', self.Ngm
+        print 'N TriggeroffsetMean: ', self.Ntom
+        
+# -----------------------------------------------------------------------------
+class fnames( object ):
+    """ organize file names of a FACT data run
+
+    """
+    
+    def __init__(self, specifier = ['012', '023', '2011', '11', '24'],
+                 rpath = '/scratch_nfs/res/bsl/',
+                 zipped = True):
+        """
+        specifier : list of strings defined as:
+            [ 'DRS calibration file', 'Data file', 'YYYY', 'MM', 'DD']
+            
+        rpath     : directory path for the results; YYYYMMDD will be appended to rpath
+        zipped    : use zipped (True) or unzipped (Data) 
+
+        """
+        
+        self.specifier = specifier
+        self.rpath     = rpath
+        self.zipped    = zipped
+        
+        self.make( self.specifier, self.rpath, self.zipped )
+
+
+    def make( self, specifier, rpath, zipped ):
+        """ create (make) the filenames
+
+        names   : dictionary of filenames, tags { 'data', 'drscal', 'results' }
+        data    : name of the data file
+        drscal  : name of the drs calibration file
+        results : radikal of file name(s) for results (to be completed  by suffixes)
+        """
+
+        self.specifier = specifier
+        
+        if zipped:
+            dpath = '/data00/fact-construction/raw/'
+            ext   = '.fits.gz'
+        else:
+            dpath = '/data03/fact-construction/raw/'
+            ext   = '.fits'
+    
+        year  = specifier[2]
+        month = specifier[3]
+        day   = specifier[4]
+        
+        yyyymmdd = year + month + day
+        dfile = specifier[1]
+        cfile = specifier[0]
+
+        rpath = rpath + yyyymmdd + '/'
+        self.rpath = rpath 
+        self.names = {}
+
+        tmp = dpath + year + '/' + month + '/' + day + '/' + yyyymmdd + '_'
+        self.names['data']  =  tmp + dfile + ext
+        self.names['drscal'] = tmp + cfile + '.drs' + ext
+        self.names['results'] =  rpath + yyyymmdd + '_' + dfile + '_' + cfile 
+
+        self.data    = self.names['data']
+        self.drscal  = self.names['drscal']
+        self.results = self.names['results']
+
+    def info( self ):
+        """ print complete filenames
+
+        """
+        
+        print 'file names:'
+        print 'data:    ', self.names['data']
+        print 'drs-cal: ', self.names['drscal']
+        print 'results: ', self.names['results']
+
+# end of class definition: fnames( object )
+class RawDataFeeder( object ):
+    """ class for wrapping multiple RawData objects
+        capable of iterating over multiple files contents as if it was one file
+    """
+    def __init__(self, filelist):
+        # actually filelist, should be a list, containing datafilepath and calibfilepath
+        # so lets check it
+        if type(filelist) != type(list()):
+            raise TypeError
+        for entry in filelist:
+            if type(entry) != type(tuple()) or len(entry) != 2:
+                raise TypeError
+    
+        self.raw_data = RawData( 
+    
+    
+class Analysis( objectg ):
+    """ basic analysis object
+        wraps a RawDataFeeder as well as 
+        input and output dictionaries
+        can write itself to disk
+    """
+    def __init__(self, outpath):
+        """ set up the analysis:
+                * write source script to self.src
+                * check if outpath is writable
+        """
+        # check if outpath is writable
+        # if not, throw IOError
+        
+        # write source
+        # I hope this works, and sys.argv[0] always contains the right filename
+        self.src = open(sys.argv[0], 'r').read()
+        pass
+    
+    def input(self, inputfile):
+        self.filelist = parse_Filelistfile(inputfile)
+        self.feeder = RawDataFeeder(self.filelist)
+        
+    def parse_Filelistfile(self, path):
+        #check inputfile type
+        if type(path) != type(str()):
+            # throw Exception TypeError or ValueError...
+        # check if file is readable, else throe IOError
+        # maybe checking if readable is enough, since os.path.methods check their input.
+        pass
+
+def _test_iter():
+    """ test for function __iter__ """
+
+#    data_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_111.fits.gz'
+#    calib_file_name =     data_file_name = 
+    data_file_name =  '/home/luster/win7/FACT/data/raw/20120114/20120114_028.fits.gz'
+    calib_file_name = '/home/luster/win7/FACT/data/raw/20120114/20120114_022.drs.fits.gz'
+    run = RawData( data_file_name, calib_file_name )
+
+    for data, scell, tt in run:
+        print 'data[0,0] = ', data[0,0], "start_cell[0] =", scell[0], "trigger type = ", tt
+
+
+if __name__ == '__main__':
+    """ tests  """
+
+    _test_iter()