source: fact/tools/pyscripts/examples/hillas_test.py@ 13227

#!/usr/bin/python -tt
#
# Example
#    * looping over RawData class object
# 
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.patches as patches

from pyfact   import RawData
from plotters import Plotter
from plotters import CamPlotter



from drs_spikes import DRSSpikes
from fir_filter import SlidingAverage
from extractor import GlobalMaxFinder
from cleaners import AmplitudeCleaner
from image_extractors import HillasParameter
from image_extractors import SimpleArea
from image_extractors import SimpleSize
import sys

data_filename = 'data/20120223_210.fits.gz'
calib_filename = 'data/20120223_206.drs.fits.gz'

run = RawData(data_filename, calib_filename, return_dict = True)
despike = DRSSpikes()
sa = SlidingAverage(8)
gmf = GlobalMaxFinder(30,230)
clean = AmplitudeCleaner(45,18)
clean.return_bool_mask = False
hillas = HillasParameter()

#p = CamPlotter('cleaned')
#p2 = CamPlotter('not')


for event in run:
    if event['trigger_type'].value == 4:
    #if True:
        print event['event_id']
        data = event['acal_data']
        unspiked_data = despike(data)
        data = sa(data)
        amp, time = gmf(data)
        survivors= clean(amp)
        
        # if nothing survived the cleaning, just go on
        if len(clean.islands) == 0:
            continue
        
        # play with the if statements here, to look at only those events you would like to analyse
        #if num_islands >= 2 and len(survivors) > 30:
        #if num_islands == 1 and len(survivors) > 20:
        if True:

            print 'calling HillasParameter'
            hillaspar = hillas(survivors, amp)
            for k in hillaspar.keys():
                if k[0] != '_':
                    print k, hillaspar[k]
                    
            #print
            #print 'delta:', hillaspar['delta']/np.pi * 180
            #print 'COG:', hillaspar['cog_euc']
            #print 'Mxx/size:', hillaspar['Mxx']/hillaspar['size']
            #print 'Mxy/size:', hillaspar['Mxy']/hillaspar['size']
            #print 'Myy/size:', hillaspar['Myy']/hillaspar['size']
            
            
            #p ( amp, survivors )
            #plt.figure( p.fig_id )
            
            # paint arroy to COG
            #plt.gca().add_patch( 
            #    patches.Polygon( 
            #    [   
            #        (hillaspar['source_pos'][0], hillaspar['source_pos'][1]), 
            #        (hillaspar['cog_euc'][0], hillaspar['cog_euc'][1]) 
            #    ] ) )
            # paint copy of x-axis through COG
            #plt.axhline( hillaspar['cog_euc'][1] )
            
            #plt.gca().add_patch( 
            #    patches.Ellipse(
            #        ( hillaspar['cog_euc'][0], hillaspar['cog_euc'][1] ),
            #        2*hillaspar['length'],
            #        2*hillaspar['width'],
            #        hillaspar['delta'] /np.pi * 180, 
            #        facecolor='none' ) )
                    
            #plt.gca().add_patch( 
            #    patches.Circle( 
            #        (hillaspar['source_pos'][0], hillaspar['source_pos'][1] ), 
            #        0.5 ) )
            
            
            
            
            #p2 (amp)
            
            #sys.exit(0)
            
            answer = raw_input('hit <Enter> to go on .... hit "q" to quit')
            if 'q' in answer:
                break
print 'good bye'