Changeset 13298 for fact/tools


Ignore:
Timestamp:
04/03/12 11:46:51 (13 years ago)
Author:
neise
Message:
made a real example out of it ... should be comprehendable and working now
File:
1 edited

Legend:

Unmodified
Added
Removed
  • fact/tools/pyscripts/examples/hillas_test.py

    r13290 r13298  
    44#    * looping over RawData class object
    55#
    6 #import matplotlib.pyplot as plt
     6import matplotlib.pyplot as plt
    77import numpy as np
    8 #import matplotlib.patches as patches
     8import matplotlib.patches as patches
    99
    1010from pyfact   import RawData
    1111#from plotters import Plotter
    12 #from plotters import CamPlotter
     12from plotters import CamPlotter
    1313
    1414
     
    3232clean = AmplitudeCleaner(45,18)
    3333clean.return_bool_mask = False
    34 hillas = HillasParameter()
     34hillas = HillasParameter( source_pos=(1,-1) )
    3535
    36 #p = CamPlotter('cleaned')
    37 #p2 = CamPlotter('not')
     36p = CamPlotter('cleaned')
    3837
    3938
     39print 'Source Position for HillasParameter class is:'
     40print hillas.source_pos
     41print 'change it in line 34 if you like'
     42
    4043for event in run:
    41     #if event['trigger_type'].value == 4:
    42     if True:
    43         print event['event_id']
     44    # trigger_type 4 means 'physical data'
     45    if event['trigger_type'].value == 4:
     46       
     47        print 'event number:', int(event['event_id'].value)
     48       
     49        #remark:
     50        # currently one needs a quite long expression for getting the event id
     51        # out of the event-dictionary:
     52        # it looks like this:         int(event['event_id'].value
     53        #
     54        # I will change it, such that it looks like this:
     55        # event['event_id']
     56        # or maybe:
     57        # event['id']
     58        #
     59        # what do you think, dear reader?? :-)
     60       
    4461        data = event['acal_data']
    4562        unspiked_data = despike(data)
     
    4865        survivors= clean(amp)
    4966       
     67        # the cleaner also sorts all pixels into groups of islands
     68        # clean.islands is a list of lists of pixel_ids
     69        # the number of lists is simply the number of islands...
     70        num_islands = len(clean.islands)
     71       
    5072        # if nothing survived the cleaning, just go on
    51         if len(clean.islands) == 0:
     73        if num_islands == 0:
    5274            continue
    5375       
    54         # play with the if statements here, to look at only those events you would like to analyse
    55         #if num_islands >= 2 and len(survivors) > 30:
    56         #if num_islands == 1 and len(survivors) > 20:
    57         if True:
    58 
    59             print 'calling HillasParameter'
     76       
     77        if num_islands == 1:
     78       
    6079            hillaspar = hillas(survivors, amp)
    61             for k in hillaspar.keys():
    62                 if k[0] != '_':
    63                     print k, hillaspar[k]
     80           
     81            print 'size', hillaspar['size']
     82            print 'alpha', hillaspar['alpha']
     83            print 'delta', hillaspar['delta']
     84            print 'width', hillaspar['width']
     85            print 'length', hillaspar['length']
     86            print 'COG in euclidic coordinates:'
     87            print hillaspar['cog_euc']
     88                       
     89           
     90            p ( amp, survivors )
     91            plt.figure( p.fig_id )
     92           
     93            # paint arrow to COG
     94            plt.gca().add_patch(
     95                patches.Polygon(
     96                [   
     97                    (hillaspar['source_pos'][0], hillaspar['source_pos'][1]),
     98                    (hillaspar['cog_euc'][0], hillaspar['cog_euc'][1])
     99                ] ) )
     100            # paint copy of x-axis through COG
     101            plt.axhline( hillaspar['cog_euc'][1] )
     102           
     103            plt.gca().add_patch(
     104                patches.Ellipse(
     105                    ( hillaspar['cog_euc'][0], hillaspar['cog_euc'][1] ),
     106                    2*hillaspar['length'],
     107                    2*hillaspar['width'],
     108                    hillaspar['delta'] /np.pi * 180,
     109                    facecolor='none' ) )
    64110                   
    65             #print
    66             #print 'delta:', hillaspar['delta']/np.pi * 180
    67             #print 'COG:', hillaspar['cog_euc']
    68             #print 'Mxx/size:', hillaspar['Mxx']/hillaspar['size']
    69             #print 'Mxy/size:', hillaspar['Mxy']/hillaspar['size']
    70             #print 'Myy/size:', hillaspar['Myy']/hillaspar['size']
     111            plt.gca().add_patch(
     112                patches.Circle(
     113                    (hillaspar['source_pos'][0], hillaspar['source_pos'][1] ),
     114                    0.5 ) )
    71115           
    72            
    73             #p ( amp, survivors )
    74             #plt.figure( p.fig_id )
    75            
    76             # paint arroy to COG
    77             #plt.gca().add_patch(
    78             #    patches.Polygon(
    79             #    [   
    80             #        (hillaspar['source_pos'][0], hillaspar['source_pos'][1]),
    81             #        (hillaspar['cog_euc'][0], hillaspar['cog_euc'][1])
    82             #    ] ) )
    83             # paint copy of x-axis through COG
    84             #plt.axhline( hillaspar['cog_euc'][1] )
    85            
    86             #plt.gca().add_patch(
    87             #    patches.Ellipse(
    88             #        ( hillaspar['cog_euc'][0], hillaspar['cog_euc'][1] ),
    89             #        2*hillaspar['length'],
    90             #        2*hillaspar['width'],
    91             #        hillaspar['delta'] /np.pi * 180,
    92             #        facecolor='none' ) )
    93                    
    94             #plt.gca().add_patch(
    95             #    patches.Circle(
    96             #        (hillaspar['source_pos'][0], hillaspar['source_pos'][1] ),
    97             #        0.5 ) )
    98            
    99            
    100            
    101            
    102             #p2 (amp)
    103            
    104             #sys.exit(0)
    105            
    106             #answer = raw_input('hit <Enter> to go on .... hit "q" to quit')
    107             #if 'q' in answer:
    108             #    break
     116            answer = raw_input('hit <Enter> to go on .... hit "q" to quit')
     117            if 'q' in answer:
     118                break
    109119print 'good bye'
Note: See TracChangeset for help on using the changeset viewer.