#!/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 to go on .... hit "q" to quit') if 'q' in answer: break print 'good bye'