Changeset 13182 for fact

Timestamp:

03/22/12 21:19:29 (13 years ago)

Author:

neise

Message:

included calculation of some more hillas parameters, still no unit tests ... no clear structure, but there is an example script in examples/hillas_test.py

File:

• fact/tools/pyscripts/pyfact/image_extractors.py (modified) (4 diffs)

fact/tools/pyscripts/pyfact/image_extractors.py

@@ -5 +5 @@
 #
 from coor import Coordinator
+import numpy as np
 from euclid import Vector2
+import sys
 
 class SimpleArea(object):
@@ -42 +44 @@
     def __init__(self):
         self.coordinator = Coordinator()
-        self.chid2coor = self.coordinator.chid2vec
+        self.chid2coor = self.coordinator.chid2coor_np
         pass
 
@@ -51 +53 @@
         self._CalculateSize()
         self._CalculateCOG()
-        
+        self._CaluculateCovarianceMatrix()
         
     def _CalculateSize(self):
@@ -69 +71 @@
         amplitude = self.amplitude
         chid2coor = self.chid2coor
+        center = self.coordinator.center
+        ex = np.array(self.coordinator.ex)
+        ey = np.array(self.coordinator.ey)
         
-        
-        cog = Vector2( 0, 0)
+        cog = np.zeros(2)
         for chid in survivors:
             cog += chid2coor[chid] * amplitude[chid]
         cog /= self.size
         self.cog = cog
-        return cog
-    def _CaluculateM(self):
-        """ M is a matrix I didn't understand yet.
-            Maybe some kind of covariance matrix...
+        print 'debug-COG=',cog, 'in hex coordinates'
+        cog_euc = cog[0] * ex + cog[1] * ey + center
+        self.cog_euc = cog_euc
+        print 'debug-COG=', cog_euc, 'in euclidic coordinates'
+        return cog, cog_euc
+        
+        
+    def _CaluculateCovarianceMatrix(self):
+        """ calculates the covariance matrix
         """
+        # shortcuts
+        survivors = self.survivors
+        amplitude = self.amplitude
+        chid2coor = self.chid2coor
+        cog_euc = self.cog_euc
+        S = self.size
+        center = self.coordinator.center
+        ex = np.array(self.coordinator.ex)
+        ey = np.array(self.coordinator.ey)
+
+        # make new 2x2 covarianz matrix
+        cov = np.zeros( (2,2) )
+        for pixel in survivors:
+            # hex coordinates
+            coor = chid2coor[pixel]
+            # make euclidic coordinates
+            coor = coor[0] * ex + coor[1] * ey + center
+            rel_coor = coor - cog_euc
+            
+            # make new 2x2 matrix
+            mi = np.ones( (2,2) )
+            # ugly ! ... 
+            mi[0][0] = rel_coor[0] * rel_coor[0]
+            mi[0][1] = rel_coor[0] * rel_coor[1]
+            mi[1][0] = rel_coor[1] * rel_coor[0]
+            mi[1][1] = rel_coor[1] * rel_coor[1]
+            mi *= amplitude[pixel]
+            #print mi , amplitude[pixel]
+            cov += mi
+        self.cov = cov
+        
+        #more funny shortcuts
+        Mxx = cov[0][0]
+        Myy = cov[1][1]
+        Mxy = cov[0][1]
+        
+        #print 'corariance matrix'
+        #print cov
+        # use the covariance matrix to calulate the number TB also cals in his thesis
+        radicand = (Myy-Mxx)**2 + 4*Mxy**2
+        k = np.sqrt(radicand) + (Myy - Mxx)
+        print 'k:', k
+        tan_delta = k / (2*Mxy )
+        if tan_delta > 0:
+            delta_in_rad = np.arctan ( tan_delta )
+        else:
+            delta_in_rad = np.arctan ( tan_delta ) + np.pi
+        delta_in_deg = delta_in_rad / np.pi * 180
+        print 'delta_in_rad' , delta_in_rad
+        print 'delta_in_deg' , delta_in_deg
+        
+        # width
+        numerator = Mxx*tan_delta**2 - k + Myy
+        denominiator = (tan_delta**2 + 1) * S
+        width = np.sqrt(numerator / denominiator)
+        
+        #length
+        numerator = Myy*tan_delta**2 + k + Mxx
+        denominiator = (tan_delta**2 + 1) * S
+        length = np.sqrt(numerator / denominiator)
+        
+        # make strange rotation matrix .. called capital R in TBs diss
+        # shortcut
+        d = delta_in_rad
+        R = np.array( ((np.cos(d), np.sin(d)),(-np.sin(d), np.cos(d))) )
+        
+        # calculate vector v between COG of shower and source position
+        # currently I don't know the cource position, so I use the center of the camera
+        v = cog_euc - np.array( [0,10] )
+        # and calculate some strange vertor r = R v 
+        r = (R*v).sum(axis=1)
+        print 'r',r
+        
+        #dist
+        dist = np.sqrt( v[0]**2 + v[1]**2 )
+        print 'dist', dist
+        
+        #alpha
+        alpha = np.arcsin( r[1] / dist)
+        alpha_in_deg = alpha / np.pi * 180
+        print 'alpha', alpha
+        print 'alpha_in_deg ', alpha_in_deg 
+        
         
 if __name__ == '__main__':