source: fact/tools/pyscripts/pyfact/image_extractors.py@ 13546

#!/usr/bin/python -tt
#
# Dominik Neise
# TU Dortmund
#
from coor import Coordinator
import numpy as np
from euclid import Vector2
import sys

class SimpleArea(object):
    """ Calculate Hillas Area in simple way.
        Sum up the number of pixels, which survived the cleaning.
        The unit is [Number of Pixel]
    """
    def __init__(self):
        # I don't know, if classes, which need no initialization
        # do still need an empty __init__ func....
        pass

    def __call__(self, survivors):
        return len(survivors)

class SimpleSize(object):
    """ Calculate Hillas Size in a very simple way
        Sum up the 'amplitude'-like value, for each survivor
    """
    def __init__(self):
        # I don't know, if classes, which need no initialization
        # do still need an empty __init__ func....
        pass

    def __call__(self, survivors, amplitude):
        size = 0
        for pixel in survivors:
            size += amplitude[pixel]
        return size

class HillasParameter(object):
    """ Calculate Hillas Parameters
        http://magic.mppmu.mpg.de/publications/theses/TBretz.pdf   p.42
        
    """
    def __init__(self, source_pos = (0,0)):
        self.coordinator = Coordinator()
        self._chid2coor = self.coordinator.chid2coor_np
        self.source_pos = source_pos

    def __call__(self, survivors, amplitude):
        self.survivors = survivors
        self.amplitude = amplitude
        
        self._CalculateSize()
        self._CalculateCOG()
        self._CalculateHillas()
        
        return self.__dict__
        
    def _CalculateSize(self):
        # shortcuts
        survivors = self.survivors
        amplitude = self.amplitude
        size = 0
        for pixel in survivors:
            size += amplitude[pixel]
            
        self.size = size
        return size
    
    def _CalculateCOG(self):
        # shortcuts
        survivors = self.survivors
        amplitude = self.amplitude
        chid2coor = self._chid2coor
        center = self.coordinator.center
        ex = np.array(self.coordinator.ex)
        ey = np.array(self.coordinator.ey)
        
        cog = np.zeros(2)
        for chid in survivors:
            cog += chid2coor[chid] * amplitude[chid]
        cog /= self.size
        self.cog = cog
        #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 _CalculateHillas(self):
        """ calculates Hillas parameters as explained in TBs dissertation on page: 42
        """
        #  define a few 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)
        source_pos = self.source_pos

        # calculate matrix elements
        
        Mxx = 0.
        Myy = 0.
        Mxy = 0.
        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
            
            
            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)
        self.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
        self.delta = delta_in_rad
        
        # width
        numerator = Mxx*tan_delta**2 - k + Myy
        denominiator = (tan_delta**2 + 1) * S
        width = np.sqrt(numerator / denominiator)
        self.width = width
        
        #length
        numerator = Myy*tan_delta**2 + k + Mxx
        denominiator = (tan_delta**2 + 1) * S
        length = np.sqrt(numerator / denominiator)
        self.length = length
        
        # 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))) )
        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 - source_pos
        # and calculate some strange vertor r = R v 
        r = (R*v).sum(axis=1)
        #print 'r',r
        self.v = v
        self.r = r
        #dist
        dist = np.sqrt( v[0]**2 + v[1]**2 )
        self.dist = dist
        #print 'dist', dist
        
        #alpha
        alpha = np.arcsin( r[1] / dist)
        alpha_in_deg = alpha / np.pi * 180
        self.alpha = alpha
        #print 'alpha', alpha
        #print 'alpha_in_deg ', alpha_in_deg 
        
        area = np.pi * width * length
        self.area = area
        
if __name__ == '__main__':
    """ no tests yet """
    print 'no test implemented yet'