| 1 | #!/usr/bin/python -itt
|
|---|
| 2 |
|
|---|
| 3 | import struct
|
|---|
| 4 | import sys
|
|---|
| 5 | import numpy as np
|
|---|
| 6 | from pprint import pprint
|
|---|
| 7 | import rlcompleter
|
|---|
| 8 | import readline
|
|---|
| 9 | readline.parse_and_bind('tab: complete')
|
|---|
| 10 | from ROOT import *
|
|---|
| 11 | import readcorsika
|
|---|
| 12 | import matplotlib.pyplot as plt
|
|---|
| 13 |
|
|---|
| 14 |
|
|---|
| 15 | def length( v ):
|
|---|
| 16 | return np.sqrt(np.dot(x,x))
|
|---|
| 17 |
|
|---|
| 18 | def matrix_times_vector( m , v):
|
|---|
| 19 | vout = v.copy()
|
|---|
| 20 | for index,line in enumerate(m):
|
|---|
| 21 | vout[index] = np.dot(line,v)
|
|---|
| 22 |
|
|---|
| 23 | return vout
|
|---|
| 24 |
|
|---|
| 25 | def make_rotation_matrix( nn, a ):
|
|---|
| 26 | R = np.zeros( (3,3) )
|
|---|
| 27 | R[0,0] = nn[0]*nn[0] * (1-np.cos(a)) + np.cos(a)
|
|---|
| 28 | R[1,0] = nn[0]*nn[1] * (1-np.cos(a)) + nn[2]*np.sin(a)
|
|---|
| 29 | R[2,0] = nn[0]*nn[2] * (1-np.cos(a)) - nn[1]*np.sin(a)
|
|---|
| 30 |
|
|---|
| 31 | R[0,1] = nn[0]*nn[1] * (1-np.cos(a)) - nn[2]*np.sin(a)
|
|---|
| 32 | R[1,1] = nn[1]*nn[1] * (1-np.cos(a)) + np.cos(a)
|
|---|
| 33 | R[2,1] = nn[1]*nn[2] * (1-np.cos(a)) + nn[0]*np.sin(a)
|
|---|
| 34 |
|
|---|
| 35 | R[0,2] = nn[0]*nn[2] * (1-np.cos(a)) + nn[1]*np.sin(a)
|
|---|
| 36 | R[1,2] = nn[1]*nn[2] * (1-np.cos(a)) - nn[0]*np.sin(a)
|
|---|
| 37 | R[2,2] = nn[2]*nn[2] * (1-np.cos(a)) + np.cos(a)
|
|---|
| 38 |
|
|---|
| 39 | return R
|
|---|
| 40 |
|
|---|
| 41 | class Thing( object ):
|
|---|
| 42 | """ Thing is just a container for the postion and the direction
|
|---|
| 43 | of something.
|
|---|
| 44 | A Thing can be a particle, or photon or something like that
|
|---|
| 45 | Or it can be a plane, like a mirror-plane or a focal-plane.
|
|---|
| 46 | Then the *dir* vector is the normal vector of the plane, and
|
|---|
| 47 | *pos* is one (possibly important) point inside the plane.
|
|---|
| 48 | """
|
|---|
| 49 | def __init__(self, pos, dir):
|
|---|
| 50 | self.pos = pos
|
|---|
| 51 | self.dir = dir
|
|---|
| 52 |
|
|---|
| 53 | def turn(self, axis, angle):
|
|---|
| 54 | """ axis might not be normalized
|
|---|
| 55 | and angle might be in degree
|
|---|
| 56 | """
|
|---|
| 57 | if length(axis) != 1.:
|
|---|
| 58 | axis /= length(axis)
|
|---|
| 59 |
|
|---|
| 60 | angle = angle / 180. *np.pi
|
|---|
| 61 |
|
|---|
| 62 |
|
|---|
| 63 |
|
|---|
| 64 | R = make_rotation_matrix( turning_axis, angle )
|
|---|
| 65 | self.pos = matrix_times_vector( R, self.pos)
|
|---|
| 66 | self.dir = matrix_times_vector( R, self.dir)
|
|---|
| 67 |
|
|---|
| 68 |
|
|---|
| 69 | def _old_turn( self, theta, phi):
|
|---|
| 70 | theta = theta/180.*np.pi
|
|---|
| 71 | phi = phi/180.*np.pi
|
|---|
| 72 |
|
|---|
| 73 | x= self.pos[0]
|
|---|
| 74 | y= self.pos[1]
|
|---|
| 75 | z= self.pos[2]
|
|---|
| 76 |
|
|---|
| 77 | vx= self.dir[0]
|
|---|
| 78 | vy= self.dir[1]
|
|---|
| 79 | vz= self.dir[2]
|
|---|
| 80 |
|
|---|
| 81 | #print vx,vy,vz
|
|---|
| 82 |
|
|---|
| 83 | #transform into spehrical coordinates
|
|---|
| 84 | print x,y,z
|
|---|
| 85 | r = length(self.pos)
|
|---|
| 86 | p = np.arctan2(y,x)
|
|---|
| 87 | t = np.arccos(z/r)
|
|---|
| 88 | print r,t,p
|
|---|
| 89 |
|
|---|
| 90 | v = length(self.dir)
|
|---|
| 91 | vp = np.arctan2(vy,vx)
|
|---|
| 92 | vt = np.arccos(vz/v)
|
|---|
| 93 |
|
|---|
| 94 | #print v,vt,vp
|
|---|
| 95 |
|
|---|
| 96 | # actual turning takes place
|
|---|
| 97 | t += theta
|
|---|
| 98 | p += phi
|
|---|
| 99 |
|
|---|
| 100 | vt += theta
|
|---|
| 101 | vp += phi
|
|---|
| 102 |
|
|---|
| 103 | #print v,vt,vp
|
|---|
| 104 |
|
|---|
| 105 | #back transform
|
|---|
| 106 |
|
|---|
| 107 | x = r * np.sin(t) * np.cos(p)
|
|---|
| 108 | y = r * np.sin(t) * np.sin(p)
|
|---|
| 109 | z = r * np.cos(t)
|
|---|
| 110 |
|
|---|
| 111 | vx = v * np.sin(vt) * np.cos(vp)
|
|---|
| 112 | vy = v * np.sin(vt) * np.sin(vp)
|
|---|
| 113 | vz = v * np.cos(vt)
|
|---|
| 114 |
|
|---|
| 115 | #print vx,vy,vz
|
|---|
| 116 |
|
|---|
| 117 | # set internal vars
|
|---|
| 118 | self.pos = np.array([x,y,z])
|
|---|
| 119 | self.dir = np.array([vx,vy,vz])
|
|---|
| 120 |
|
|---|
| 121 |
|
|---|
| 122 | def __repr__( self ):
|
|---|
| 123 | return "%s(%r)" % (self.__class__, self.__dict__)
|
|---|
| 124 |
|
|---|
| 125 | class Mirror( Thing ):
|
|---|
| 126 | def __init__(self, index, pos, normal_vector, focal_length, hex_size ):
|
|---|
| 127 | super(Mirror,self).__init__(pos, normal_vector)
|
|---|
| 128 | self.index = index
|
|---|
| 129 | self.focal_length = focal_length
|
|---|
| 130 | self.hex_size = hex_size
|
|---|
| 131 |
|
|---|
| 132 | def __repr__( self ):
|
|---|
| 133 | return "%s(%r)" % (self.__class__, self.__dict__)
|
|---|
| 134 |
|
|---|
| 135 |
|
|---|
| 136 |
|
|---|
| 137 |
|
|---|
| 138 |
|
|---|
| 139 |
|
|---|
| 140 | def read_reflector_definition_file( filename ):
|
|---|
| 141 | """
|
|---|
| 142 | """
|
|---|
| 143 | mirrors = []
|
|---|
| 144 |
|
|---|
| 145 | f = open( filename )
|
|---|
| 146 | for index, line in enumerate(f):
|
|---|
| 147 | if line[0] == '#':
|
|---|
| 148 | continue
|
|---|
| 149 | line = line.split()
|
|---|
| 150 | if len(line) < 8:
|
|---|
| 151 | continue
|
|---|
| 152 | #print line
|
|---|
| 153 |
|
|---|
| 154 | # first 3 colums in the file are x,y,z coordinates of the center
|
|---|
| 155 | # of this mirror in cm, I guess
|
|---|
| 156 | pos = np.array(map(float, line[0:3]))
|
|---|
| 157 |
|
|---|
| 158 | # the next 3 elements are the elements of the normal vector
|
|---|
| 159 | # should be normalized already, so the unit is of no importance.
|
|---|
| 160 | normal_vector = np.array(map(float,line[3:6]))
|
|---|
| 161 |
|
|---|
| 162 | # focal length of this mirror in mm
|
|---|
| 163 | focal_length = float(line[6])
|
|---|
| 164 |
|
|---|
| 165 | # size of the hexagonal shaped facette mirror, measured as the radius
|
|---|
| 166 | # of the hexagons *inner* circle.
|
|---|
| 167 | hex_size = float(line[8])
|
|---|
| 168 | mirror = Mirror( index, pos, normal_vector, focal_length, hex_size )
|
|---|
| 169 |
|
|---|
| 170 | mirrors.append(mirror)
|
|---|
| 171 |
|
|---|
| 172 | return mirrors
|
|---|
| 173 |
|
|---|
| 174 |
|
|---|
| 175 | def reflect_photon( photon, mirrors):
|
|---|
| 176 | """ finds out:
|
|---|
| 177 | which mirror is hit by photon
|
|---|
| 178 | and where
|
|---|
| 179 | and in which angle relative to mirror
|
|---|
| 180 | """
|
|---|
| 181 |
|
|---|
| 182 | # the line defined by the photon is used to find the intersection point
|
|---|
| 183 | # with the plane of each facette mirror. Then I check,
|
|---|
| 184 | # if the intersection point lies within the limits of the facette mirrors
|
|---|
| 185 | # hexagonal boundaries.
|
|---|
| 186 | # If this is the case I have found the mirror, which is hit, and
|
|---|
| 187 | # can calculate:
|
|---|
| 188 | # the distance of the intersection point from the center of the facette mirror
|
|---|
| 189 | # and the angle relative to the mirror (normal or plane not sure yet)
|
|---|
| 190 |
|
|---|
| 191 | for mirror in mirrors:
|
|---|
| 192 | #facette mirror plane, defined as n . x = d1 . n
|
|---|
| 193 | n = mirror.dir
|
|---|
| 194 | d1 = mirror.pos
|
|---|
| 195 |
|
|---|
| 196 | # line of photon defined as r = lambda * v + d2
|
|---|
| 197 | v = photon.dir
|
|---|
| 198 | d2 = photon.pos
|
|---|
| 199 |
|
|---|
| 200 | # the intersection coordinates are found by solving
|
|---|
| 201 | # n . (lambda * v + d2) - d1 . n == 0, for lambda=lambda_0
|
|---|
| 202 | # and then the intersection is: i = lambda_0 * v + d2
|
|---|
| 203 | #
|
|---|
| 204 | # putting int in another form:
|
|---|
| 205 | # solve:
|
|---|
| 206 | # lambda * n.v + n.d2 - n.d1 == 0
|
|---|
| 207 | # or
|
|---|
| 208 | # lambda_0 = n.(d1-d2) / n.v
|
|---|
| 209 |
|
|---|
| 210 | # FIXME: if one of the two dot-products is very small,
|
|---|
| 211 | # we shuold take special care maybe
|
|---|
| 212 | # if n.(d1-d2) is very small, this means that the starting point of
|
|---|
| 213 | # the photon is already nearly in the plane, so lambda_0 is expected to
|
|---|
| 214 | # be very small ... erm .. maybe this is actually not a special case
|
|---|
| 215 | # but very good.
|
|---|
| 216 | # of n.v is very small, this means the patch of the photon is nearly
|
|---|
| 217 | # parallel to the plane, so the error ob lambda_0 might be very large.
|
|---|
| 218 | # in addition, this might just tell us, that the mirror is hit under
|
|---|
| 219 | # strange circumstances ... so its not a good candidate and we can already go on.
|
|---|
| 220 | lambda_0 = (np.dot(n,(d1-d2)) / np.dot(n,v))
|
|---|
| 221 |
|
|---|
| 222 | #intersection between line and plane
|
|---|
| 223 | i = lambda_0 * v + d2
|
|---|
| 224 |
|
|---|
| 225 | # I want the distance beween i and d1 so I can already from the distance find
|
|---|
| 226 | # out if this is our candidate.
|
|---|
| 227 | distance = np.sqrt(((i-d1)*(i-d1)).sum())
|
|---|
| 228 |
|
|---|
| 229 | #print "photon pos:", d2, "\t dir:", v/length(v)
|
|---|
| 230 | #print "mirror pos:", d1, "\t dir:", n/length(n)
|
|---|
| 231 |
|
|---|
| 232 | #print "lambda_0", lambda_0
|
|---|
| 233 | #print "intersection :", i
|
|---|
| 234 | #print "distance:",distance
|
|---|
| 235 |
|
|---|
| 236 | if distance <= mirror.hex_size/2.:
|
|---|
| 237 | break
|
|---|
| 238 | else:
|
|---|
| 239 | mirror = None
|
|---|
| 240 |
|
|---|
| 241 | if not mirror is None:
|
|---|
| 242 | photon.mirror_index = mirror.index
|
|---|
| 243 | photon.mirror_intersection = i
|
|---|
| 244 | photon.mirror_center_distance = distance
|
|---|
| 245 | #print "mirror found:", mirror.index ,
|
|---|
| 246 | #print "distance", distance
|
|---|
| 247 | # now I have to find out, if the photon is not only in the
|
|---|
| 248 | # right distance but actually has hit the mirror.
|
|---|
| 249 | # this I do like this
|
|---|
| 250 | # i-d1 is a vector in the mirror plane pointing from d1 to the intersection point i.
|
|---|
| 251 | # if I know turn the entire mirror plane so it lies withing the x-y-plane
|
|---|
| 252 | # by applying a simple turning-matrix, then each vector inside the plane with turn into
|
|---|
| 253 | # a nice x,x vector.
|
|---|
| 254 | # now I assume, that the hexagon is "pointing" lets say to into y direction
|
|---|
| 255 | # so I can e.g. say:
|
|---|
| 256 | # x has to be between -30.3 and +30.3 and y has to be
|
|---|
| 257 | # between 35 - m * |x| and -35 + m * |x| ... pretty simple.
|
|---|
| 258 | # maybe one can leave the turning aside, but I like that I can imagine it very nicely
|
|---|
| 259 | #
|
|---|
| 260 | #
|
|---|
| 261 | # I don't do this yet .. since I don't know by heart how a turning matrix looks :-)
|
|---|
| 262 | # so I just simulate round mirrors
|
|---|
| 263 | ######################################################################
|
|---|
| 264 |
|
|---|
| 265 |
|
|---|
| 266 | # next step, since I know the intersection point, is the new direction.
|
|---|
| 267 | # So I need the normal of the mirror in the intersection point.
|
|---|
| 268 | # since the normal of every mirror is alway pointing to the camera center
|
|---|
| 269 | # this is not difficult.
|
|---|
| 270 |
|
|---|
| 271 | normal_at_intersection = (mirror_alignmen_point.pos - i) / length(mirror_alignmen_point.pos - i)
|
|---|
| 272 | #print "normal_at_intersection",normal_at_intersection
|
|---|
| 273 |
|
|---|
| 274 | angle = np.arccos(np.dot( v, normal_at_intersection) / (length(v) * length(normal_at_intersection)))
|
|---|
| 275 | photon.angle_to_mirror_normal = angle
|
|---|
| 276 | #print "angle:", angle/np.pi*180., "deg"
|
|---|
| 277 |
|
|---|
| 278 |
|
|---|
| 279 | # okay, now I have the intersection *i*,
|
|---|
| 280 | # the old direction of the photon *v*
|
|---|
| 281 | # and the normalvector at the intersection.
|
|---|
| 282 | ######################################################################
|
|---|
| 283 | ######################################################################
|
|---|
| 284 | # I do this now differently.
|
|---|
| 285 | # I will mirror the "point" at the tip of *v* at the line created by
|
|---|
| 286 | # the normalvector at the intersection and the intersection.
|
|---|
| 287 | # this will gibe me a mirrored_point *mp* and the vector from *i* to *mp*
|
|---|
| 288 | # is the *new_direction* it should even be normalized.
|
|---|
| 289 |
|
|---|
| 290 | # 1. step: create plane through the "tip" of *v* and the normal_at_intersection.
|
|---|
| 291 | # 2. step: find crossingpoint on line through *i* and the normal_at_intersection,
|
|---|
| 292 | # 3. step: vector from "tip" of *v* to crossingpoint times 2 points to
|
|---|
| 293 | # the "tip" of *mirrored_v*
|
|---|
| 294 |
|
|---|
| 295 | # plane: n_plane_3 . r = p_plane_3 . n_plane_3
|
|---|
| 296 | # p_plane_3 = i+v
|
|---|
| 297 | # n_plane_3 = normal_at_intersection
|
|---|
| 298 |
|
|---|
| 299 | # line: r = lambda_3 * v_line_3 + p_line_3
|
|---|
| 300 | # p_line_3 = i
|
|---|
| 301 | # v_line_3 = normal_at_intersection
|
|---|
| 302 |
|
|---|
| 303 | # create crossing: n_plane_3 . (lambda_3 * v_line_3 + p_line_3) = p_plane_3 . n_plane_3
|
|---|
| 304 | # <=> lambda_3 = (p_plane_3 - p_line_3 ).n_plane_3 / n_plane_3 . v_line_3
|
|---|
| 305 | # <=> lambda_3 = (i+v - i).normal_at_intersection / normal_at_intersection . normal_at_intersection
|
|---|
| 306 | # <=> lambda_3 = v.normal_at_intersection
|
|---|
| 307 |
|
|---|
| 308 | lambda_3 = np.dot(v, normal_at_intersection)
|
|---|
| 309 | #print "lambda_3", lambda_3
|
|---|
| 310 | crossing_point_3 = lambda_3 * normal_at_intersection + i
|
|---|
| 311 | #print "crossing_point_3", crossing_point_3
|
|---|
| 312 |
|
|---|
| 313 | from_tip_of_v_to_crossing_point_3 = crossing_point_3 - (i+v)
|
|---|
| 314 |
|
|---|
| 315 | tip_of_mirrored_v = i+v+ 2*from_tip_of_v_to_crossing_point_3
|
|---|
| 316 |
|
|---|
| 317 | new_direction = tip_of_mirrored_v - i
|
|---|
| 318 |
|
|---|
| 319 | #print "new_direction",new_direction
|
|---|
| 320 | #print "old direction", v
|
|---|
| 321 | photon.new_direction = new_direction
|
|---|
| 322 | ######################################################################
|
|---|
| 323 | ######################################################################
|
|---|
| 324 | """
|
|---|
| 325 | # both directions form a plane, and when I turn the old *v* by
|
|---|
| 326 | # twice the angle between *v* and *normal_at_intersection*
|
|---|
| 327 | # inside this plane then I get the new direction of the photon.
|
|---|
| 328 |
|
|---|
| 329 | # so lets first get the normal of the reflection plane
|
|---|
| 330 | normal_of_reflection_plane =np.cross( v, normal_at_intersection)
|
|---|
| 331 |
|
|---|
| 332 | print length(normal_of_reflection_plane), "should be one"
|
|---|
| 333 | print length(v), "should be one"
|
|---|
| 334 | print length(normal_at_intersection), "should be one"
|
|---|
| 335 | print np.dot(v, normal_at_intersection), "should *NOT* be zero"
|
|---|
| 336 | print np.dot(v, normal_of_reflection_plane), "should be zero"
|
|---|
| 337 | print np.dot(normal_at_intersection, normal_of_reflection_plane), "should be zero"
|
|---|
| 338 |
|
|---|
| 339 | angle = np.arccos(np.dot( v, normal_at_intersection) / (length(v) * length(normal_at_intersection)))
|
|---|
| 340 | photon.angle_to_mirror_normal = angle
|
|---|
| 341 | print "angle:", angle/np.pi*180., "deg"
|
|---|
| 342 |
|
|---|
| 343 | # the rotation matrix for the rotation of *v* around normal_of_reflection_plane is
|
|---|
| 344 | R = make_rotation_matrix( normal_of_reflection_plane, 2*angle )
|
|---|
| 345 |
|
|---|
| 346 | print "R"
|
|---|
| 347 | pprint(R)
|
|---|
| 348 |
|
|---|
| 349 | new_direction = matrix_times_vector( R, v)
|
|---|
| 350 | photon.new_direction = new_direction
|
|---|
| 351 |
|
|---|
| 352 | print "old direction", v, length(v)
|
|---|
| 353 | print "new direction", new_direction, length(new_direction)
|
|---|
| 354 | print "mirror center", mirror.pos
|
|---|
| 355 | print "interception point", i
|
|---|
| 356 | print "center of focal plane", focal_plane.pos
|
|---|
| 357 | """
|
|---|
| 358 |
|
|---|
| 359 | # new the photon has a new direction *new_direction* and is starting
|
|---|
| 360 | # from the intersection point *i*
|
|---|
| 361 | # now I want to find out where there focal plane is hit.
|
|---|
| 362 | # So I have to repeat the stuff from up there
|
|---|
| 363 |
|
|---|
| 364 | #print "np.dot(focal_plane.dir,new_direction))", np.dot(focal_plane.dir,new_direction)
|
|---|
| 365 |
|
|---|
| 366 | lambda_1 = (np.dot(focal_plane.dir ,(focal_plane.pos - i)) / np.dot(focal_plane.dir,new_direction))
|
|---|
| 367 |
|
|---|
| 368 | #print "lambda_1", lambda_1
|
|---|
| 369 | focal_plane_spot = lambda_1 * new_direction + i
|
|---|
| 370 | #print "focal_plane_spot",focal_plane_spot
|
|---|
| 371 | photon.hit = True
|
|---|
| 372 | photon.focal_plane_pos = focal_plane_spot - focal_plane.pos
|
|---|
| 373 |
|
|---|
| 374 | #print "distance from focal plane center=", length(focal_plane_spot-focal_plane.pos)
|
|---|
| 375 | else:
|
|---|
| 376 | photon.hit = False
|
|---|
| 377 |
|
|---|
| 378 |
|
|---|
| 379 | return photon
|
|---|
| 380 |
|
|---|
| 381 |
|
|---|
| 382 |
|
|---|
| 383 | def reflect_photon_old( photon, mirrors):
|
|---|
| 384 | """ finds out:
|
|---|
| 385 | which mirror is hit by photon
|
|---|
| 386 | and where
|
|---|
| 387 | and in which angle relative to mirror
|
|---|
| 388 | """
|
|---|
| 389 |
|
|---|
| 390 | # the line defined by the photon is used to find the intersection point
|
|---|
| 391 | # with the plane of each facette mirror. Then I check,
|
|---|
| 392 | # if the intersection point lies within the limits of the facette mirrors
|
|---|
| 393 | # hexagonal boundaries.
|
|---|
| 394 | # If this is the case I have found the mirror, which is hit, and
|
|---|
| 395 | # can calculate:
|
|---|
| 396 | # the distance of the intersection point from the center of the facette mirror
|
|---|
| 397 | # and the angle relative to the mirror (normal or plane not sure yet)
|
|---|
| 398 |
|
|---|
| 399 | for mirror in mirrors:
|
|---|
| 400 | #facette mirror plane, defined as n . x = d1 . n
|
|---|
| 401 | n = mirror.dir
|
|---|
| 402 | d1 = mirror.pos
|
|---|
| 403 |
|
|---|
| 404 | # line of photon defined as r = lambda * v + d2
|
|---|
| 405 | v = photon.dir
|
|---|
| 406 | d2 = photon.pos
|
|---|
| 407 |
|
|---|
| 408 | # the intersection coordinates are found by solving
|
|---|
| 409 | # n . (lambda * v + d2) - d1 . n == 0, for lambda=lambda_0
|
|---|
| 410 | # and then the intersection is: i = lambda_0 * v + d2
|
|---|
| 411 | #
|
|---|
| 412 | # putting int in another form:
|
|---|
| 413 | # solve:
|
|---|
| 414 | # lambda * n.v + n.d2 - n.d1 == 0
|
|---|
| 415 | # or
|
|---|
| 416 | # lambda_0 = n.(d1-d2) / n.v
|
|---|
| 417 |
|
|---|
| 418 | # FIXME: if one of the two dot-products is very small,
|
|---|
| 419 | # we shuold take special care maybe
|
|---|
| 420 | # if n.(d1-d2) is very small, this means that the starting point of
|
|---|
| 421 | # the photon is already nearly in the plane, so lambda_0 is expected to
|
|---|
| 422 | # be very small ... erm .. maybe this is actually not a special case
|
|---|
| 423 | # but very good.
|
|---|
| 424 | # of n.v is very small, this means the patch of the photon is nearly
|
|---|
| 425 | # parallel to the plane, so the error ob lambda_0 might be very large.
|
|---|
| 426 | # in addition, this might just tell us, that the mirror is hit under
|
|---|
| 427 | # strange circumstances ... so its not a good candidate and we can already go on.
|
|---|
| 428 | lambda_0 = (np.dot(n,(d1-d2)) / np.dot(n,v))
|
|---|
| 429 |
|
|---|
| 430 | #intersection between line and plane
|
|---|
| 431 | i = lambda_0 * v + d2
|
|---|
| 432 |
|
|---|
| 433 | # I want the distance beween i and d1 so I can already from the distance find
|
|---|
| 434 | # out if this is our candidate.
|
|---|
| 435 | distance = np.sqrt(((i-d1)*(i-d1)).sum())
|
|---|
| 436 |
|
|---|
| 437 | #print "photon pos:", d2, "\t dir:", v/length(v)
|
|---|
| 438 | #print "mirror pos:", d1, "\t dir:", n/length(n)
|
|---|
| 439 |
|
|---|
| 440 | #print "lambda_0", lambda_0
|
|---|
| 441 | #print "intersection :", i
|
|---|
| 442 | #print "distance:",distance
|
|---|
| 443 |
|
|---|
| 444 | if distance <= mirror.hex_size/2.:
|
|---|
| 445 | break
|
|---|
| 446 | else:
|
|---|
| 447 | mirror = None
|
|---|
| 448 |
|
|---|
| 449 | if not mirror is None:
|
|---|
| 450 | photon.mirror_index = mirror.index
|
|---|
| 451 | photon.mirror_intersection = i
|
|---|
| 452 | photon.mirror_center_distance = distance
|
|---|
| 453 | #print "mirror found:", mirror.index ,
|
|---|
| 454 | #print "distance", distance
|
|---|
| 455 | # now I have to find out, if the photon is not only in the
|
|---|
| 456 | # right distance but actually has hit the mirror.
|
|---|
| 457 | # this I do like this
|
|---|
| 458 | # i-d1 is a vector in the mirror plane pointing from d1 to the intersection point i.
|
|---|
| 459 | # if I know turn the entire mirror plane so it lies withing the x-y-plane
|
|---|
| 460 | # by applying a simple turning-matrix, then each vector inside the plane with turn into
|
|---|
| 461 | # a nice x,x vector.
|
|---|
| 462 | # now I assume, that the hexagon is "pointing" lets say to into y direction
|
|---|
| 463 | # so I can e.g. say:
|
|---|
| 464 | # x has to be between -30.3 and +30.3 and y has to be
|
|---|
| 465 | # between 35 - m * |x| and -35 + m * |x| ... pretty simple.
|
|---|
| 466 | # maybe one can leave the turning aside, but I like that I can imagine it very nicely
|
|---|
| 467 | #
|
|---|
| 468 | #
|
|---|
| 469 | # I don't do this yet .. since I don't know by heart how a turning matrix looks :-)
|
|---|
| 470 | # so I just simulate round mirrors
|
|---|
| 471 | ######################################################################
|
|---|
| 472 |
|
|---|
| 473 |
|
|---|
| 474 | # next step, since I know the intersection point, is the new direction.
|
|---|
| 475 | # So I need the normal of the mirror in the intersection point.
|
|---|
| 476 | # since the normal of every mirror is alway pointing to the camera center
|
|---|
| 477 | # this is not difficult.
|
|---|
| 478 |
|
|---|
| 479 | normal_at_intersection = (mirror_alignmen_point.pos - i) / length(mirror_alignmen_point.pos - i)
|
|---|
| 480 |
|
|---|
| 481 | # okay, now I have the intersection *i*,
|
|---|
| 482 | # the old direction of the photon *v*
|
|---|
| 483 | # and the normalvector at the intersection.
|
|---|
| 484 | # both directions form a plane, and when I turn the old *v* by
|
|---|
| 485 | # twice the angle between *v* and *normal_at_intersection*
|
|---|
| 486 | # inside this plane then I get the new direction of the photon.
|
|---|
| 487 |
|
|---|
| 488 | # so lets first get the normal of the reflection plane
|
|---|
| 489 | normal_of_reflection_plane =np.cross( v, normal_at_intersection)
|
|---|
| 490 |
|
|---|
| 491 | #print length(normal_of_reflection_plane)
|
|---|
| 492 | #print length(v)
|
|---|
| 493 | #print length(normal_at_intersection)
|
|---|
| 494 | #print np.dot(v, normal_at_intersection)
|
|---|
| 495 | #print np.dot(v, normal_of_reflection_plane)
|
|---|
| 496 | #print np.dot(normal_at_intersection, normal_of_reflection_plane)
|
|---|
| 497 |
|
|---|
| 498 | angle = np.arccos(np.dot( v, normal_at_intersection) / (length(v) * length(normal_at_intersection)))
|
|---|
| 499 | photon.angle_to_mirror_normal = angle
|
|---|
| 500 | #print "angle:", angle/np.pi*180., "deg"
|
|---|
| 501 |
|
|---|
| 502 | # the rotation matrix for the rotation of *v* around normal_of_reflection_plane is
|
|---|
| 503 | R = make_rotation_matrix( normal_of_reflection_plane, 2*angle )
|
|---|
| 504 |
|
|---|
| 505 |
|
|---|
| 506 | new_direction = matrix_times_vector( R, v)
|
|---|
| 507 | photon.new_direction = new_direction
|
|---|
| 508 |
|
|---|
| 509 | #print "old direction", v, length(v)
|
|---|
| 510 | #print "new direction", new_direction, length(new_direction)
|
|---|
| 511 | #print "mirror center", mirror.pos
|
|---|
| 512 | #print "interception point", i
|
|---|
| 513 | #print "center of focal plane", focal_plane.pos
|
|---|
| 514 |
|
|---|
| 515 |
|
|---|
| 516 | # new the photon has a new direction *new_direction* and is starting
|
|---|
| 517 | # from the intersection point *i*
|
|---|
| 518 | # now I want to find out where there focal plane is hit.
|
|---|
| 519 | # So I have to repeat the stuff from up there
|
|---|
| 520 |
|
|---|
| 521 | lambda_1 = (np.dot(focal_plane.dir ,(focal_plane.pos - i)) / np.dot(focal_plane.dir,new_direction))
|
|---|
| 522 | #print "lambda_1", lambda_1
|
|---|
| 523 | focal_plane_spot = lambda_1 * new_direction + i
|
|---|
| 524 | photon.hit = True
|
|---|
| 525 | photon.focal_plane_pos = focal_plane_spot - focal_plane.pos
|
|---|
| 526 |
|
|---|
| 527 | #print "distance from focal plane center=", length(focal_plane_spot-focal_plane.pos)
|
|---|
| 528 | else:
|
|---|
| 529 | photon.hit = False
|
|---|
| 530 |
|
|---|
| 531 |
|
|---|
| 532 | return photon
|
|---|
| 533 |
|
|---|
| 534 | class Photon( Thing ):
|
|---|
| 535 | """ a photon has not only the direction and position, which a Thing has.
|
|---|
| 536 | but it also has a wavelength and a "time" and a "mother_particle_id"
|
|---|
| 537 |
|
|---|
| 538 | the photon constructor understands the 10-element 1D-np.array
|
|---|
| 539 | which is stored inside a run.event.photons 2D-np.array
|
|---|
| 540 | """
|
|---|
| 541 | def __init__(self, photon_definition_array ):
|
|---|
| 542 | """ the *photon_definition_array* pda contains:
|
|---|
| 543 | pda[0] - encoded info
|
|---|
| 544 | pda[1:3] - x,y position in cm
|
|---|
| 545 | pda[3:5] - u,v cosines to x,y axis --> so called direction cosines
|
|---|
| 546 | pda[5] - time since first interaction [ns]
|
|---|
| 547 | pda[6] - height of production in cm
|
|---|
| 548 | pda[7] - j ??
|
|---|
| 549 | pda[8] - imov ??
|
|---|
| 550 | pda[9] - wavelength [nm]
|
|---|
| 551 | """
|
|---|
| 552 | pda = photon_definition_array
|
|---|
| 553 | pos = np.array([pda[1],pda[2],0.])
|
|---|
| 554 | dir = np.array([pda[3],pda[4], np.sqrt(1-pda[3]**2-pda[4]**2) ])
|
|---|
| 555 | super(Photon,self).__init__(pos, dir)
|
|---|
| 556 |
|
|---|
| 557 | self.wavelength = pda[9]
|
|---|
| 558 | self.time = pda[5]
|
|---|
| 559 |
|
|---|
| 560 | def __repr__( self ):
|
|---|
| 561 | return "%s(%r)" % (self.__class__, self.__dict__)
|
|---|
| 562 |
|
|---|
| 563 |
|
|---|
| 564 | if __name__ == '__main__':
|
|---|
| 565 | mirrors = read_reflector_definition_file( "030/fact-reflector.txt" )
|
|---|
| 566 |
|
|---|
| 567 |
|
|---|
| 568 | focal_plane = Thing( pos=np.array([0.,0.,978.132/2.]), # center of focal_plane
|
|---|
| 569 | dir=np.array([0., 0., 1.]) ) # direction of view
|
|---|
| 570 | focal_plane.size = 20 # diameter in cm
|
|---|
| 571 |
|
|---|
| 572 |
|
|---|
| 573 | mirror_alignmen_point = Thing( pos=np.array([0.,0.,978.132]), # center of focal_plane
|
|---|
| 574 | dir=np.array([0., 0., 1.]) ) # direction of view
|
|---|
| 575 |
|
|---|
| 576 |
|
|---|
| 577 |
|
|---|
| 578 | # turn the telescope
|
|---|
| 579 | turning_axis = np.array([1,0,0])
|
|---|
| 580 | angle = -3
|
|---|
| 581 |
|
|---|
| 582 | for mirror in mirrors:
|
|---|
| 583 | mirror.turn( turning_axis, angle)
|
|---|
| 584 |
|
|---|
| 585 | focal_plane.turn( turning_axis, angle)
|
|---|
| 586 | mirror_alignmen_point.turn( turning_axis, angle)
|
|---|
| 587 |
|
|---|
| 588 |
|
|---|
| 589 | turning_axis = np.array([0,1,0])
|
|---|
| 590 | angle = 30
|
|---|
| 591 |
|
|---|
| 592 | for mirror in mirrors:
|
|---|
| 593 | mirror.turn( turning_axis, angle)
|
|---|
| 594 |
|
|---|
| 595 | focal_plane.turn( turning_axis, angle)
|
|---|
| 596 | mirror_alignmen_point.turn( turning_axis, angle)
|
|---|
| 597 |
|
|---|
| 598 |
|
|---|
| 599 |
|
|---|
| 600 |
|
|---|
| 601 | run = readcorsika.read_corsika_file("cer")
|
|---|
| 602 |
|
|---|
| 603 | li=[]
|
|---|
| 604 | event_counter = 0
|
|---|
| 605 | for event in run.events:
|
|---|
| 606 | print event_counter
|
|---|
| 607 | event_counter += 1
|
|---|
| 608 | event.photons_who_hit = []
|
|---|
| 609 | for photon in event.photons:
|
|---|
| 610 | # photon is a 1D-nd.array with 10 elements
|
|---|
| 611 | # make Photon instance from this horrible array
|
|---|
| 612 | photon[1:3] -= np.array(event.info_dict['core_loc'])
|
|---|
| 613 | photon = Photon(photon)
|
|---|
| 614 |
|
|---|
| 615 | photon = reflect_photon( photon, mirrors )
|
|---|
| 616 |
|
|---|
| 617 | if photon.hit:
|
|---|
| 618 | li.append( photon.angle_to_mirror_normal / np.pi * 180.)
|
|---|
| 619 | event.photons_who_hit.append(photon)
|
|---|
| 620 |
|
|---|
| 621 | if event_counter > 100:
|
|---|
| 622 | break
|
|---|
| 623 |
|
|---|
| 624 | #plt.ion()
|
|---|
| 625 | #fig1 = plt.figure()
|
|---|
| 626 | #fig2 = plt.figure()
|
|---|
| 627 |
|
|---|
| 628 | #plt.hist( np.array(li) , bins=100)
|
|---|
| 629 | plt.hold(False)
|
|---|
| 630 |
|
|---|
| 631 |
|
|---|
| 632 | g = TGraph2D()
|
|---|
| 633 | h = TH2F("h","title",120,-60,60,120,-60,60)
|
|---|
| 634 | c = 0
|
|---|
| 635 | for ev in run.events:
|
|---|
| 636 |
|
|---|
| 637 | ground_dirs = []
|
|---|
| 638 | focal_positions = []
|
|---|
| 639 | if not hasattr(ev, "photons_who_hit"):
|
|---|
| 640 | continue
|
|---|
| 641 | for ph in ev.photons_who_hit:
|
|---|
| 642 | #mi = ph.mirror_intersection
|
|---|
| 643 | new = Thing( pos=ph.focal_plane_pos, dir = ph.new_direction)
|
|---|
| 644 | new.turn( np.array([0,1,0]), -30 )
|
|---|
| 645 | mi = new.pos
|
|---|
| 646 | #print mi
|
|---|
| 647 | h.Fill(mi[0], mi[1])
|
|---|
| 648 | g.SetPoint( c, mi[0], mi[1], mi[2])
|
|---|
| 649 | c += 1
|
|---|
| 650 | ground_dirs.append(ph.dir)
|
|---|
| 651 | focal_positions.append(new.pos)
|
|---|
| 652 |
|
|---|
| 653 | ground_dirs = np.array(ground_dirs)
|
|---|
| 654 | focal_positions = np.array(focal_positions)
|
|---|
| 655 | #print ground_dirs.shape, focal_positions.shape
|
|---|
| 656 | #plt.figure(fig1.number)
|
|---|
| 657 | #plt.plot( ground_dirs[:,0], ground_dirs[:,1], '.')
|
|---|
| 658 | #plt.title("ground directions")
|
|---|
| 659 | #plt.figure(fig2.number)
|
|---|
| 660 | #plt.plot( focal_positions[:,0], focal_positions[:,1], '.')
|
|---|
| 661 | #plt.title("focal positions")
|
|---|
| 662 |
|
|---|
| 663 |
|
|---|
| 664 | #raw_input()
|
|---|
| 665 |
|
|---|
| 666 | g.SetMarkerStyle(20)
|
|---|
| 667 | #g.Draw("pcol")
|
|---|
| 668 | h.Draw("colz")
|
|---|
| 669 |
|
|---|
| 670 |
|
|---|
| 671 |
|
|---|
| 672 |
|
|---|
| 673 | """
|
|---|
| 674 | for m in mirrors:
|
|---|
| 675 | mi = m.pos
|
|---|
| 676 | #mi = ph.focal_plane_pos
|
|---|
| 677 | g.SetPoint( c, mi[0], mi[1], mi[2])
|
|---|
| 678 | c += 1
|
|---|
| 679 | """
|
|---|
| 680 |
|
|---|
| 681 |
|
|---|
| 682 |
|
|---|
| 683 |
|
|---|
| 684 |
|
|---|
| 685 |
|
|---|
| 686 |
|
|---|
