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

#!/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( source_pos=(1,-1) )

p = CamPlotter('cleaned')


print 'Source Position for HillasParameter class is:'
print hillas.source_pos
print 'change it in line 34 if you like'

for event in run:
    # trigger_type 4 means 'physical data'
    if event['trigger_type'].value == 4:
        
        print 'event number:', int(event['event_id'].value)
        
        #remark:
        # currently one needs a quite long expression for getting the event id
        # out of the event-dictionary:
        # it looks like this:         int(event['event_id'].value
        #
        # I will change it, such that it looks like this:
        # event['event_id']
        # or maybe:
        # event['id']
        #
        # what do you think, dear reader?? :-)
        
        data = event['acal_data']
        unspiked_data = despike(data)
        data = sa(data)
        amp, time = gmf(data)
        survivors= clean(amp)
        
        # the cleaner also sorts all pixels into groups of islands
        # clean.islands is a list of lists of pixel_ids 
        # the number of lists is simply the number of islands...
        num_islands = len(clean.islands)
        
        # if nothing survived the cleaning, just go on
        if num_islands == 0:
            continue
        
        
        if num_islands == 1:
        
            hillaspar = hillas(survivors, amp)
            
            print 'size', hillaspar['size']
            print 'alpha', hillaspar['alpha']
            print 'delta', hillaspar['delta']
            print 'width', hillaspar['width']
            print 'length', hillaspar['length']
            print 'COG in euclidic coordinates:'
            print hillaspar['cog_euc']
                        
            
            p ( amp, survivors )
            plt.figure( p.fig_id )
            
            # paint arrow 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 ) )
            
            answer = raw_input('hit <Enter> to go on .... hit "q" to quit')
            if 'q' in answer:
                break
print 'good bye'