Changeset 13298

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:

• fact/tools/pyscripts/examples/hillas_test.py (modified) (3 diffs)

fact/tools/pyscripts/examples/hillas_test.py

@@ -4 +4 @@
 #    * looping over RawData class object
 # 
-#import matplotlib.pyplot as plt
+import matplotlib.pyplot as plt
 import numpy as np
-#import matplotlib.patches as patches
+import matplotlib.patches as patches
 
 from pyfact   import RawData
 #from plotters import Plotter
-#from plotters import CamPlotter
+from plotters import CamPlotter
 
 
@@ -32 +32 @@
 clean = AmplitudeCleaner(45,18)
 clean.return_bool_mask = False
-hillas = HillasParameter()
+hillas = HillasParameter( source_pos=(1,-1) )
 
-#p = CamPlotter('cleaned')
-#p2 = CamPlotter('not')
+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:
-    #if event['trigger_type'].value == 4:
-    if True:
-        print event['event_id']
+    # 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)
@@ -48 +65 @@
         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 len(clean.islands) == 0:
+        if num_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'
+        
+        if num_islands == 1:
+        
             hillaspar = hillas(survivors, amp)
-            for k in hillaspar.keys():
-                if k[0] != '_':
-                    print k, hillaspar[k]
+            
+            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' ) )
                     
-            #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']
+            plt.gca().add_patch( 
+                patches.Circle( 
+                    (hillaspar['source_pos'][0], hillaspar['source_pos'][1] ), 
+                    0.5 ) )
             
-            
-            #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
+            answer = raw_input('hit <Enter> to go on .... hit "q" to quit')
+            if 'q' in answer:
+                break
 print 'good bye'