source: fact/tools/pyscripts/sandbox/vogler/LP_analysis_from_Patrick__example.py

#!/usr/bin/python -tti
#
# Lightpulser Analysis Test Script
# Written by DN but
# but based on an idea from PV 
# This script is intended for testing purposes, 
# and has no real purpose so far....
#
print 'importing libraries ... takes a sec .. sorry'

# 2012 04 17
#calibfilename = '/fact/raw/2012/04/17/20120417_003.drs.fits.gz'    # NROI 300, pedestal
calibfilename = '/fact/raw/2012/04/17/20120417_033.drs.fits.gz'    # NROI 300, pedestal

#datafilename = '/fact/raw/2012/04/17/20120417_015.fits.gz'         # NROI 300, LP_ext
#datafilename = '/fact/raw/2012/04/17/20120417_021.fits.gz'         # NROI 300, LP_ext
#datafilename = '/fact/raw/2012/04/17/20120417_036.fits.gz'         # NROI 300, LP_ext
datafilename = '/fact/raw/2012/04/17/20120417_042.fits.gz'         # NROI 300, LP_ext

# 2012 01 25
#datafilename = '/fact/raw/2012/01/25/20120125_042.fits.gz'     # light-pulser-ext
#  = '/fact/raw/2012/01/25/20120125_095.fits.gz'     # pedestal
#  = '/fact/raw/2012/01/25/20120125_094.fits.gz'     # pedestal
#  = '/fact/raw/2012/01/25/20120125_093.fits.gz'     # pedestal
#datafilename = '/fact/raw/2012/01/25/20120125_075.fits.gz'
#calibfilename = '/fact/raw/2012/01/25/20120125_088.drs.fits.gz'


from pyfact import RawData
import numpy as np
import sys
from plotters import Plotter
from plotters import CamPlotter

print 'creating RawData object...'
run = RawData(datafilename, calibfilename, return_dict=True)

print "... looping..."
print ' over ', run.nevents, 'events  ... please wait and stay calm :-)'

# create variables as 'None' so inside the loop the variables can be created 
# with the right shape
data_average = None
extracted_average = None

for event in run:
    data = event['acal_data'] # this data is already in the shape number_of_pixel x ROI

    # create an array to store the 'average event', the shape derived from the shape of the data
    if data_average == None:
        data_average = np.zeros( data.shape )
    data_average += data


    # A first test of Adrians idea of primitive signal extractor
    # signal and background
    
    # axis=1 means, do not sum over all, but only over the slices, the output is a numpy array with len=numpix
    extracted_signal = data[ : , 90:100 ].sum( axis=1 )  
    
    # background subtraction
    extracted_signal -= data[ : , 15:25 ].sum( axis=1 )

    if extracted_average == None:
        extracted_average = np.zeros( extracted_signal.shape )
    extracted_average += extracted_signal

print "Looped ... "

# After the loop we should divide 
#   the data_average and 
#   the extracted_average
# by the number of events
data_average /= run.nevents
extracted_average /= 10 * run.nevents  # the 10 is hardcoded here ... this is bad coding ... but it works :-)

# I guess data_average_run was intended to contain the mean over all pixel
# the name is a bit strange ...
# axis=0 means, to calculate the mean over all pixel.
# so outcome is a 1D array of length = ROI
data_average_run = data_average.mean( axis=0 )

myplotter = Plotter('this is myplotter', xlabel='time in slices', ylabel='amplitude calibrated data ... in mV')
mycamplotter = CamPlotter('this is mycamplotter',  vmax=350)  # for Lightpulser data
#mycamplotter = CamPlotter('titel of the plot',  map_file_path = '../map_dn.txt', vmin=-1.2, vmax=1.2)    # for pedestal data

myplotter( data_average[22], 'pix 22, average' )
mycamplotter( extracted_average )

print 
print 'the script was started in "interactive mode" by using the cmdline option -tti in the first line of the script'
print 'you have still all the variables from the script available here and can go on testing interactively....'
print 
print 'eg type:'
print " p = Plotter('test')"
print " p(data_average_run) "
print 
print ' the RawData instance called run, was not yet deleted, since in interactive mode this makes no sense ... you might still need it :-)'