| 1 | #!/usr/bin/python -tt
|
|---|
| 2 | #
|
|---|
| 3 | # Dominik Neise
|
|---|
| 4 | # TU Dortmund
|
|---|
| 5 | #
|
|---|
| 6 | from coor import Coordinator
|
|---|
| 7 | import numpy as np
|
|---|
| 8 | from euclid import Vector2
|
|---|
| 9 | import sys
|
|---|
| 10 |
|
|---|
| 11 | class SimpleArea(object):
|
|---|
| 12 | """ Calculate Hillas Area in simple way.
|
|---|
| 13 | Sum up the number of pixels, which survived the cleaning.
|
|---|
| 14 | The unit is [Number of Pixel]
|
|---|
| 15 | """
|
|---|
| 16 | def __init__(self):
|
|---|
| 17 | # I don't know, if classes, which need no initialization
|
|---|
| 18 | # do still need an empty __init__ func....
|
|---|
| 19 | pass
|
|---|
| 20 |
|
|---|
| 21 | def __call__(self, survivors):
|
|---|
| 22 | return len(survivors)
|
|---|
| 23 |
|
|---|
| 24 | class SimpleSize(object):
|
|---|
| 25 | """ Calculate Hillas Size in a very simple way
|
|---|
| 26 | Sum up the 'amplitude'-like value, for each survivor
|
|---|
| 27 | """
|
|---|
| 28 | def __init__(self):
|
|---|
| 29 | # I don't know, if classes, which need no initialization
|
|---|
| 30 | # do still need an empty __init__ func....
|
|---|
| 31 | pass
|
|---|
| 32 |
|
|---|
| 33 | def __call__(self, survivors, amplitude):
|
|---|
| 34 | size = 0
|
|---|
| 35 | for pixel in survivors:
|
|---|
| 36 | size += amplitude[pixel]
|
|---|
| 37 | return size
|
|---|
| 38 |
|
|---|
| 39 | class HillasParameter(object):
|
|---|
| 40 | """ Calculate Hillas Parameters
|
|---|
| 41 | http://magic.mppmu.mpg.de/publications/theses/TBretz.pdf p.42
|
|---|
| 42 |
|
|---|
| 43 | """
|
|---|
| 44 | def __init__(self, source_pos = (0,0)):
|
|---|
| 45 | self.coordinator = Coordinator()
|
|---|
| 46 | self._chid2coor = self.coordinator.chid2coor_np
|
|---|
| 47 | self.source_pos = source_pos
|
|---|
| 48 |
|
|---|
| 49 | def __call__(self, survivors, amplitude):
|
|---|
| 50 | self.survivors = survivors
|
|---|
| 51 | self.amplitude = amplitude
|
|---|
| 52 |
|
|---|
| 53 | self._CalculateSize()
|
|---|
| 54 | self._CalculateCOG()
|
|---|
| 55 | self._CalculateHillas()
|
|---|
| 56 |
|
|---|
| 57 | return self.__dict__
|
|---|
| 58 |
|
|---|
| 59 | def _CalculateSize(self):
|
|---|
| 60 | # shortcuts
|
|---|
| 61 | survivors = self.survivors
|
|---|
| 62 | amplitude = self.amplitude
|
|---|
| 63 | size = 0
|
|---|
| 64 | for pixel in survivors:
|
|---|
| 65 | size += amplitude[pixel]
|
|---|
| 66 |
|
|---|
| 67 | self.size = size
|
|---|
| 68 | return size
|
|---|
| 69 |
|
|---|
| 70 | def _CalculateCOG(self):
|
|---|
| 71 | # shortcuts
|
|---|
| 72 | survivors = self.survivors
|
|---|
| 73 | amplitude = self.amplitude
|
|---|
| 74 | chid2coor = self._chid2coor
|
|---|
| 75 | center = self.coordinator.center
|
|---|
| 76 | ex = np.array(self.coordinator.ex)
|
|---|
| 77 | ey = np.array(self.coordinator.ey)
|
|---|
| 78 |
|
|---|
| 79 | cog = np.zeros(2)
|
|---|
| 80 | for chid in survivors:
|
|---|
| 81 | cog += chid2coor[chid] * amplitude[chid]
|
|---|
| 82 | cog /= self.size
|
|---|
| 83 | self.cog = cog
|
|---|
| 84 | #print 'debug-COG=',cog, 'in hex coordinates'
|
|---|
| 85 | cog_euc = cog[0] * ex + cog[1] * ey + center
|
|---|
| 86 | self.cog_euc = cog_euc
|
|---|
| 87 | #print 'debug-COG=', cog_euc, 'in euclidic coordinates'
|
|---|
| 88 | return cog, cog_euc
|
|---|
| 89 |
|
|---|
| 90 |
|
|---|
| 91 | def _CalculateHillas(self):
|
|---|
| 92 | """ calculates Hillas parameters as explained in TBs dissertation on page: 42
|
|---|
| 93 | """
|
|---|
| 94 | # define a few shortcuts
|
|---|
| 95 | survivors = self.survivors
|
|---|
| 96 | amplitude = self.amplitude
|
|---|
| 97 | chid2coor = self._chid2coor
|
|---|
| 98 | cog_euc = self.cog_euc
|
|---|
| 99 | S = self.size
|
|---|
| 100 | center = self.coordinator.center
|
|---|
| 101 | ex = np.array(self.coordinator.ex)
|
|---|
| 102 | ey = np.array(self.coordinator.ey)
|
|---|
| 103 | source_pos = self.source_pos
|
|---|
| 104 |
|
|---|
| 105 | # calculate matrix elements
|
|---|
| 106 |
|
|---|
| 107 | Mxx = 0.
|
|---|
| 108 | Myy = 0.
|
|---|
| 109 | Mxy = 0.
|
|---|
| 110 | for pixel in survivors:
|
|---|
| 111 | # hex coordinates
|
|---|
| 112 | coor = chid2coor[pixel]
|
|---|
| 113 | # make euclidic coordinates
|
|---|
| 114 | coor = coor[0] * ex + coor[1] * ey + center
|
|---|
| 115 | rel_coor = coor - cog_euc
|
|---|
| 116 |
|
|---|
| 117 |
|
|---|
| 118 | mxx = rel_coor[0] * rel_coor[0] * amplitude[pixel]
|
|---|
| 119 | mxy = rel_coor[0] * rel_coor[1] * amplitude[pixel]
|
|---|
| 120 | myy = rel_coor[1] * rel_coor[1] * amplitude[pixel]
|
|---|
| 121 |
|
|---|
| 122 | Mxx += mxx
|
|---|
| 123 | Myy += myy
|
|---|
| 124 | Mxy += mxy
|
|---|
| 125 | self.Mxx = Mxx
|
|---|
| 126 | self.Mxy = Mxy
|
|---|
| 127 | self.Myy = Myy
|
|---|
| 128 |
|
|---|
| 129 |
|
|---|
| 130 | # use the covariance matrix to calulate the number TB also cals in his thesis
|
|---|
| 131 | radicand = (Myy-Mxx)**2 + 4*Mxy**2
|
|---|
| 132 | k = np.sqrt(radicand) + (Myy - Mxx)
|
|---|
| 133 | self.k = k
|
|---|
| 134 | tan_delta = k / (2*Mxy )
|
|---|
| 135 | if tan_delta > 0:
|
|---|
| 136 | delta_in_rad = np.arctan ( tan_delta )
|
|---|
| 137 | else:
|
|---|
| 138 | delta_in_rad = np.arctan ( tan_delta ) + np.pi
|
|---|
| 139 | delta_in_deg = delta_in_rad / np.pi * 180
|
|---|
| 140 | #print 'delta_in_rad' , delta_in_rad
|
|---|
| 141 | #print 'delta_in_deg' , delta_in_deg
|
|---|
| 142 | self.delta = delta_in_rad
|
|---|
| 143 |
|
|---|
| 144 | # width
|
|---|
| 145 | numerator = Mxx*tan_delta**2 - k + Myy
|
|---|
| 146 | denominiator = (tan_delta**2 + 1) * S
|
|---|
| 147 | width = np.sqrt(numerator / denominiator)
|
|---|
| 148 | self.width = width
|
|---|
| 149 |
|
|---|
| 150 | #length
|
|---|
| 151 | numerator = Myy*tan_delta**2 + k + Mxx
|
|---|
| 152 | denominiator = (tan_delta**2 + 1) * S
|
|---|
| 153 | length = np.sqrt(numerator / denominiator)
|
|---|
| 154 | self.length = length
|
|---|
| 155 |
|
|---|
| 156 | # make strange rotation matrix .. called capital R in TBs diss
|
|---|
| 157 | # shortcut
|
|---|
| 158 | d = delta_in_rad
|
|---|
| 159 | R = np.array( ((np.cos(d), np.sin(d)),(-np.sin(d), np.cos(d))) )
|
|---|
| 160 | self.R = R
|
|---|
| 161 | # calculate vector v between COG of shower and source position
|
|---|
| 162 | # currently I don't know the cource position, so I use the center of the camera
|
|---|
| 163 |
|
|---|
| 164 | v = cog_euc - source_pos
|
|---|
| 165 | # and calculate some strange vertor r = R v
|
|---|
| 166 | r = (R*v).sum(axis=1)
|
|---|
| 167 | #print 'r',r
|
|---|
| 168 | self.v = v
|
|---|
| 169 | self.r = r
|
|---|
| 170 | #dist
|
|---|
| 171 | dist = np.sqrt( v[0]**2 + v[1]**2 )
|
|---|
| 172 | self.dist = dist
|
|---|
| 173 | #print 'dist', dist
|
|---|
| 174 |
|
|---|
| 175 | #alpha
|
|---|
| 176 | alpha = np.arcsin( r[1] / dist)
|
|---|
| 177 | alpha_in_deg = alpha / np.pi * 180
|
|---|
| 178 | self.alpha = alpha
|
|---|
| 179 | #print 'alpha', alpha
|
|---|
| 180 | #print 'alpha_in_deg ', alpha_in_deg
|
|---|
| 181 |
|
|---|
| 182 | area = np.pi * width * length
|
|---|
| 183 | self.area = area
|
|---|
| 184 |
|
|---|
| 185 | if __name__ == '__main__':
|
|---|
| 186 | """ no tests yet """
|
|---|
| 187 | print 'no test implemented yet' |
|---|
