Changeset 13229

Timestamp:

03/26/12 13:32:44 (13 years ago)

Author:

neise

Message:

still testing

File:

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

fact/tools/pyscripts/pyfact/image_extractors.py

@@ -42 +42 @@
         
     """
-    def __init__(self):
+    def __init__(self, source_pos = (0,0)):
         self.coordinator = Coordinator()
-        self.chid2coor = self.coordinator.chid2coor_np
-        pass
+        self._chid2coor = self.coordinator.chid2coor_np
+        self.source_pos = source_pos
 
     def __call__(self, survivors, amplitude):
@@ -53 +53 @@
         self._CalculateSize()
         self._CalculateCOG()
-        self._CaluculateCovarianceMatrix()
+        self._CalculateHillas()
+        
+        return self.__dict__
         
     def _CalculateSize(self):
@@ -70 +72 @@
         survivors = self.survivors
         amplitude = self.amplitude
-        chid2coor = self.chid2coor
+        chid2coor = self._chid2coor
         center = self.coordinator.center
         ex = np.array(self.coordinator.ex)
@@ -87 +89 @@
         
         
-    def _CaluculateCovarianceMatrix(self):
-        """ calculates the covariance matrix
+    def _CalculateHillas(self):
+        """ calculates Hillas parameters as explained in TBs dissertation on page: 42
         """
-        # shortcuts
+        #  define a few shortcuts
         survivors = self.survivors
         amplitude = self.amplitude
-        chid2coor = self.chid2coor
+        chid2coor = self._chid2coor
         cog_euc = self.cog_euc
         S = self.size
@@ -99 +101 @@
         ex = np.array(self.coordinator.ex)
         ey = np.array(self.coordinator.ey)
+        source_pos = self.source_pos
 
-        # make new 2x2 covarianz matrix
-        cov = np.zeros( (2,2) )
+        # calculate matrix elements
+        
+        Mxx = 0.
+        Myy = 0.
+        Mxy = 0.
         for pixel in survivors:
             # hex coordinates
@@ -109 +115 @@
             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
+            
+            mxx = rel_coor[0] * rel_coor[0] * amplitude[pixel]
+            mxy = rel_coor[0] * rel_coor[1] * amplitude[pixel]
+            myy = rel_coor[1] * rel_coor[1] * amplitude[pixel]
+
+            Mxx += mxx
+            Myy += myy
+            Mxy += mxy
+        self.Mxx = Mxx
+        self.Mxy = Mxy
+        self.Myy = Myy
+
+
         # 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
+        self.k = k
         tan_delta = k / (2*Mxy )
         if tan_delta > 0:
@@ -138 +138 @@
             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
+        #print 'delta_in_rad' , delta_in_rad
+        #print 'delta_in_deg' , delta_in_deg
+        self.delta = delta_in_rad
         
         # width
@@ -145 +146 @@
         denominiator = (tan_delta**2 + 1) * S
         width = np.sqrt(numerator / denominiator)
+        self.width = width
         
         #length
@@ -150 +152 @@
         denominiator = (tan_delta**2 + 1) * S
         length = np.sqrt(numerator / denominiator)
+        self.length = length
         
         # make strange rotation matrix .. called capital R in TBs diss
@@ -155 +158 @@
         d = delta_in_rad
         R = np.array( ((np.cos(d), np.sin(d)),(-np.sin(d), np.cos(d))) )
-        
+        self.R = R
         # 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] )
+
+        v = cog_euc - source_pos
         # and calculate some strange vertor r = R v 
         r = (R*v).sum(axis=1)
-        print 'r',r
-        
+        #print 'r',r
+        self.v = v
+        self.r = r
         #dist
         dist = np.sqrt( v[0]**2 + v[1]**2 )
-        print 'dist', dist
+        self.dist = dist
+        #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 
+        self.alpha = alpha
+        #print 'alpha', alpha
+        #print 'alpha_in_deg ', alpha_in_deg 
         
+        area = np.pi * width * length
+        self.area = area
         
 if __name__ == '__main__':