| 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 | from Turnable import *
|
|---|
| 15 |
|
|---|
| 16 | import simulation_helpers as sh
|
|---|
| 17 |
|
|---|
| 18 | class Mirror( Turnable ):
|
|---|
| 19 | """ Mirror description/abstraction class
|
|---|
| 20 |
|
|---|
| 21 | from Turnable:
|
|---|
| 22 | turn( axis, angle) - turn the mirror around *axis* by *angle*
|
|---|
| 23 |
|
|---|
| 24 | *pos* - the position of the center of the mirror plane
|
|---|
| 25 |
|
|---|
| 26 | *turnables*:
|
|---|
| 27 | *pos*
|
|---|
| 28 | *dir*
|
|---|
| 29 |
|
|---|
| 30 | """
|
|---|
| 31 | def __init__(self, index, pos, normal_vector, focal_length, hex_size ):
|
|---|
| 32 | self.pos = pos
|
|---|
| 33 | self.dir = normal_vector
|
|---|
| 34 |
|
|---|
| 35 | # list for Turnable.turn()
|
|---|
| 36 | self._turnables = [ "pos", "dir"]
|
|---|
| 37 |
|
|---|
| 38 | self.index = index
|
|---|
| 39 | self.focal_length = focal_length
|
|---|
| 40 | self.hex_size = hex_size
|
|---|
| 41 |
|
|---|
| 42 | # standard __repr__
|
|---|
| 43 | def __repr__( self ):
|
|---|
| 44 | return "%s(%r)" % (self.__class__, self.__dict__)
|
|---|
| 45 |
|
|---|
| 46 | class Photon( Turnable ):
|
|---|
| 47 | """ Photon description/abstraction class
|
|---|
| 48 |
|
|---|
| 49 | from Turnable:
|
|---|
| 50 | turn( axis, angle) - turn the mirror around *axis* by *angle*
|
|---|
| 51 |
|
|---|
| 52 | *pos* - the position of the photon in 3D
|
|---|
| 53 | *dir* - the direction of flight in 3D
|
|---|
| 54 | *time*- time since 2st interaction [ns]
|
|---|
| 55 | *wavelength* - wavelength in [nm]
|
|---|
| 56 | *mother* - particle ID of mother particle
|
|---|
| 57 |
|
|---|
| 58 | *turnables*:
|
|---|
| 59 | *pos*
|
|---|
| 60 | *dir*
|
|---|
| 61 |
|
|---|
| 62 | """
|
|---|
| 63 | def __init__(self, photon_definition_array ):
|
|---|
| 64 | """ Construct Photon form 10-element 1D-np.array
|
|---|
| 65 |
|
|---|
| 66 | the photon constructor understands the 10-element 1D-np.array
|
|---|
| 67 | which is stored inside a run.event.photons 2D-np.array
|
|---|
| 68 |
|
|---|
| 69 | the *photon_definition_array* pda contains:
|
|---|
| 70 | pda[0] - encoded info
|
|---|
| 71 | pda[1:3] - x,y position in cm
|
|---|
| 72 | pda[3:5] - u,v cosines to x,y axis --> so called direction cosines
|
|---|
| 73 | pda[5] - time since first interaction [ns]
|
|---|
| 74 | pda[6] - height of production in cm
|
|---|
| 75 | pda[7] - j ??
|
|---|
| 76 | pda[8] - imov ??
|
|---|
| 77 | pda[9] - wavelength [nm]
|
|---|
| 78 | """
|
|---|
| 79 | pda = photon_definition_array
|
|---|
| 80 |
|
|---|
| 81 | self.pos = np.array([pda[1],pda[2],0.])
|
|---|
| 82 | self.dir = np.array([pda[3],pda[4], np.sqrt(1-pda[3]**2-pda[4]**2) ])
|
|---|
| 83 | self._turnables = ("pos", "dir")
|
|---|
| 84 |
|
|---|
| 85 | self.wavelength = pda[9]
|
|---|
| 86 | self.time = pda[5]
|
|---|
| 87 |
|
|---|
| 88 | # standard __repr__
|
|---|
| 89 | def __repr__( self ):
|
|---|
| 90 | return "%s(%r)" % (self.__class__, self.__dict__)
|
|---|
| 91 |
|
|---|
| 92 |
|
|---|
| 93 | class Focal_Plane( Turnable ):
|
|---|
| 94 | """
|
|---|
| 95 | """
|
|---|
| 96 | def __init__(self, pos, dir, size ):
|
|---|
| 97 | self.pos = pos
|
|---|
| 98 | self.dir = dir
|
|---|
| 99 | self.size = size
|
|---|
| 100 | self._turnables = ["pos", "dir"]
|
|---|
| 101 |
|
|---|
| 102 | class Point( Turnable ):
|
|---|
| 103 | """
|
|---|
| 104 | """
|
|---|
| 105 | def __init__(self, pos):
|
|---|
| 106 | self.pos = pos
|
|---|
| 107 | self._turnables = ["pos"]
|
|---|
| 108 |
|
|---|
| 109 |
|
|---|
| 110 |
|
|---|
| 111 | def read_reflector_definition_file( filename ):
|
|---|
| 112 | """
|
|---|
| 113 | """
|
|---|
| 114 | mirrors = []
|
|---|
| 115 |
|
|---|
| 116 | f = open( filename )
|
|---|
| 117 | for index, line in enumerate(f):
|
|---|
| 118 | if line[0] == '#':
|
|---|
| 119 | continue
|
|---|
| 120 | line = line.split()
|
|---|
| 121 | if len(line) < 8:
|
|---|
| 122 | continue
|
|---|
| 123 | #print line
|
|---|
| 124 |
|
|---|
| 125 | # first 3 colums in the file are x,y,z coordinates of the center
|
|---|
| 126 | # of this mirror in cm, I guess
|
|---|
| 127 | pos = np.array(map(float, line[0:3]))
|
|---|
| 128 |
|
|---|
| 129 | # the next 3 elements are the elements of the normal vector
|
|---|
| 130 | # should be normalized already, so the unit is of no importance.
|
|---|
| 131 | normal_vector = np.array(map(float,line[3:6]))
|
|---|
| 132 |
|
|---|
| 133 | # focal length of this mirror in mm
|
|---|
| 134 | focal_length = float(line[6])
|
|---|
| 135 |
|
|---|
| 136 | # size of the hexagonal shaped facette mirror, measured as the radius
|
|---|
| 137 | # of the hexagons *inner* circle.
|
|---|
| 138 | hex_size = float(line[8])
|
|---|
| 139 | mirror = Mirror( index, pos, normal_vector, focal_length, hex_size )
|
|---|
| 140 |
|
|---|
| 141 | mirrors.append(mirror)
|
|---|
| 142 |
|
|---|
| 143 | return mirrors
|
|---|
| 144 |
|
|---|
| 145 |
|
|---|
| 146 |
|
|---|
| 147 |
|
|---|
| 148 | def reflect_photon( photon, mirrors):
|
|---|
| 149 | """ finds out:
|
|---|
| 150 | which mirror is hit by photon
|
|---|
| 151 | and where
|
|---|
| 152 | and in which angle relative to mirror
|
|---|
| 153 | """
|
|---|
| 154 |
|
|---|
| 155 | # the line defined by the photon is used to find the intersection point
|
|---|
| 156 | # with the plane of each facette mirror. Then I check,
|
|---|
| 157 | # if the intersection point lies within the limits of the facette mirrors
|
|---|
| 158 | # hexagonal boundaries.
|
|---|
| 159 | # If this is the case I have found the mirror, which is hit, and
|
|---|
| 160 | # can calculate:
|
|---|
| 161 | # the distance of the intersection point from the center of the facette mirror
|
|---|
| 162 | # and the angle relative to the mirror (normal or plane not sure yet)
|
|---|
| 163 |
|
|---|
| 164 | for mirror in mirrors:
|
|---|
| 165 | #facette mirror plane, defined as n . x = d1 . n
|
|---|
| 166 | n = mirror.dir
|
|---|
| 167 | d1 = mirror.pos
|
|---|
| 168 |
|
|---|
| 169 | # line of photon defined as r = lambda * v + d2
|
|---|
| 170 | v = photon.dir
|
|---|
| 171 | d2 = photon.pos
|
|---|
| 172 |
|
|---|
| 173 | # the intersection coordinates are found by solving
|
|---|
| 174 | # n . (lambda * v + d2) - d1 . n == 0, for lambda=lambda_0
|
|---|
| 175 | # and then the intersection is: i = lambda_0 * v + d2
|
|---|
| 176 | #
|
|---|
| 177 | # putting int in another form:
|
|---|
| 178 | # solve:
|
|---|
| 179 | # lambda * n.v + n.d2 - n.d1 == 0
|
|---|
| 180 | # or
|
|---|
| 181 | # lambda_0 = n.(d1-d2) / n.v
|
|---|
| 182 |
|
|---|
| 183 | # FIXME: if one of the two dot-products is very small,
|
|---|
| 184 | # we shuold take special care maybe
|
|---|
| 185 | # if n.(d1-d2) is very small, this means that the starting point of
|
|---|
| 186 | # the photon is already nearly in the plane, so lambda_0 is expected to
|
|---|
| 187 | # be very small ... erm .. maybe this is actually not a special case
|
|---|
| 188 | # but very good.
|
|---|
| 189 | # of n.v is very small, this means the patch of the photon is nearly
|
|---|
| 190 | # parallel to the plane, so the error ob lambda_0 might be very large.
|
|---|
| 191 | # in addition, this might just tell us, that the mirror is hit under
|
|---|
| 192 | # strange circumstances ... so its not a good candidate and we can already go on.
|
|---|
| 193 | lambda_0 = (np.dot(n,(d1-d2)) / np.dot(n,v))
|
|---|
| 194 |
|
|---|
| 195 | #intersection between line and plane
|
|---|
| 196 | i = lambda_0 * v + d2
|
|---|
| 197 |
|
|---|
| 198 | # I want the distance beween i and d1 so I can already from the distance find
|
|---|
| 199 | # out if this is our candidate.
|
|---|
| 200 | distance = np.sqrt(((i-d1)*(i-d1)).sum())
|
|---|
| 201 |
|
|---|
| 202 | #print "photon pos:", d2, "\t dir:", v/length(v)
|
|---|
| 203 | #print "mirror pos:", d1, "\t dir:", n/length(n)
|
|---|
| 204 |
|
|---|
| 205 | #print "lambda_0", lambda_0
|
|---|
| 206 | #print "intersection :", i
|
|---|
| 207 | #print "distance:",distance
|
|---|
| 208 |
|
|---|
| 209 | if distance <= mirror.hex_size/2.:
|
|---|
| 210 | break
|
|---|
| 211 | else:
|
|---|
| 212 | mirror = None
|
|---|
| 213 |
|
|---|
| 214 | if not mirror is None:
|
|---|
| 215 | photon.mirror_index = mirror.index
|
|---|
| 216 | photon.mirror_intersection = i
|
|---|
| 217 | photon.mirror_center_distance = distance
|
|---|
| 218 | #print "mirror found:", mirror.index ,
|
|---|
| 219 | #print "distance", distance
|
|---|
| 220 | # now I have to find out, if the photon is not only in the
|
|---|
| 221 | # right distance but actually has hit the mirror.
|
|---|
| 222 | # this I do like this
|
|---|
| 223 | # i-d1 is a vector in the mirror plane pointing from d1 to the intersection point i.
|
|---|
| 224 | # if I know turn the entire mirror plane so it lies withing the x-y-plane
|
|---|
| 225 | # by applying a simple turning-matrix, then each vector inside the plane with turn into
|
|---|
| 226 | # a nice x,x vector.
|
|---|
| 227 | # now I assume, that the hexagon is "pointing" lets say to into y direction
|
|---|
| 228 | # so I can e.g. say:
|
|---|
| 229 | # x has to be between -30.3 and +30.3 and y has to be
|
|---|
| 230 | # between 35 - m * |x| and -35 + m * |x| ... pretty simple.
|
|---|
| 231 | # maybe one can leave the turning aside, but I like that I can imagine it very nicely
|
|---|
| 232 | #
|
|---|
| 233 | #
|
|---|
| 234 | # I don't do this yet .. since I don't know by heart how a turning matrix looks :-)
|
|---|
| 235 | # so I just simulate round mirrors
|
|---|
| 236 | ######################################################################
|
|---|
| 237 |
|
|---|
| 238 |
|
|---|
| 239 | # next step, since I know the intersection point, is the new direction.
|
|---|
| 240 | # So I need the normal of the mirror in the intersection point.
|
|---|
| 241 | # since the normal of every mirror is alway pointing to the camera center
|
|---|
| 242 | # this is not difficult.
|
|---|
| 243 |
|
|---|
| 244 | normal_at_intersection = (mirror_alignmen_point.pos - i) / sh.length(mirror_alignmen_point.pos - i)
|
|---|
| 245 | #print "normal_at_intersection",normal_at_intersection
|
|---|
| 246 |
|
|---|
| 247 | angle = np.arccos(np.dot( v, normal_at_intersection) / (sh.length(v) * sh.length(normal_at_intersection)))
|
|---|
| 248 | photon.angle_to_mirror_normal = angle
|
|---|
| 249 | #print "angle:", angle/np.pi*180., "deg"
|
|---|
| 250 |
|
|---|
| 251 |
|
|---|
| 252 | # okay, now I have the intersection *i*,
|
|---|
| 253 | # the old direction of the photon *v*
|
|---|
| 254 | # and the normalvector at the intersection.
|
|---|
| 255 | ######################################################################
|
|---|
| 256 | ######################################################################
|
|---|
| 257 | # I do this now differently.
|
|---|
| 258 | # I will mirror the "point" at the tip of *v* at the line created by
|
|---|
| 259 | # the normalvector at the intersection and the intersection.
|
|---|
| 260 | # this will gibe me a mirrored_point *mp* and the vector from *i* to *mp*
|
|---|
| 261 | # is the *new_direction* it should even be normalized.
|
|---|
| 262 |
|
|---|
| 263 | # 1. step: create plane through the "tip" of *v* and the normal_at_intersection.
|
|---|
| 264 | # 2. step: find crossingpoint on line through *i* and the normal_at_intersection,
|
|---|
| 265 | # 3. step: vector from "tip" of *v* to crossingpoint times 2 points to
|
|---|
| 266 | # the "tip" of *mirrored_v*
|
|---|
| 267 |
|
|---|
| 268 | # plane: n_plane_3 . r = p_plane_3 . n_plane_3
|
|---|
| 269 | # p_plane_3 = i+v
|
|---|
| 270 | # n_plane_3 = normal_at_intersection
|
|---|
| 271 |
|
|---|
| 272 | # line: r = lambda_3 * v_line_3 + p_line_3
|
|---|
| 273 | # p_line_3 = i
|
|---|
| 274 | # v_line_3 = normal_at_intersection
|
|---|
| 275 |
|
|---|
| 276 | # create crossing: n_plane_3 . (lambda_3 * v_line_3 + p_line_3) = p_plane_3 . n_plane_3
|
|---|
| 277 | # <=> lambda_3 = (p_plane_3 - p_line_3 ).n_plane_3 / n_plane_3 . v_line_3
|
|---|
| 278 | # <=> lambda_3 = (i+v - i).normal_at_intersection / normal_at_intersection . normal_at_intersection
|
|---|
| 279 | # <=> lambda_3 = v.normal_at_intersection
|
|---|
| 280 |
|
|---|
| 281 | lambda_3 = np.dot(v, normal_at_intersection)
|
|---|
| 282 | #print "lambda_3", lambda_3
|
|---|
| 283 | crossing_point_3 = lambda_3 * normal_at_intersection + i
|
|---|
| 284 | #print "crossing_point_3", crossing_point_3
|
|---|
| 285 |
|
|---|
| 286 | from_tip_of_v_to_crossing_point_3 = crossing_point_3 - (i+v)
|
|---|
| 287 |
|
|---|
| 288 | tip_of_mirrored_v = i+v+ 2*from_tip_of_v_to_crossing_point_3
|
|---|
| 289 |
|
|---|
| 290 | new_direction = tip_of_mirrored_v - i
|
|---|
| 291 |
|
|---|
| 292 | #print "new_direction",new_direction
|
|---|
| 293 | #print "old direction", v
|
|---|
| 294 | photon.new_direction = new_direction
|
|---|
| 295 | ######################################################################
|
|---|
| 296 | ######################################################################
|
|---|
| 297 | """
|
|---|
| 298 | # both directions form a plane, and when I turn the old *v* by
|
|---|
| 299 | # twice the angle between *v* and *normal_at_intersection*
|
|---|
| 300 | # inside this plane then I get the new direction of the photon.
|
|---|
| 301 |
|
|---|
| 302 | # so lets first get the normal of the reflection plane
|
|---|
| 303 | normal_of_reflection_plane =np.cross( v, normal_at_intersection)
|
|---|
| 304 |
|
|---|
| 305 | print length(normal_of_reflection_plane), "should be one"
|
|---|
| 306 | print length(v), "should be one"
|
|---|
| 307 | print length(normal_at_intersection), "should be one"
|
|---|
| 308 | print np.dot(v, normal_at_intersection), "should *NOT* be zero"
|
|---|
| 309 | print np.dot(v, normal_of_reflection_plane), "should be zero"
|
|---|
| 310 | print np.dot(normal_at_intersection, normal_of_reflection_plane), "should be zero"
|
|---|
| 311 |
|
|---|
| 312 | angle = np.arccos(np.dot( v, normal_at_intersection) / (length(v) * length(normal_at_intersection)))
|
|---|
| 313 | photon.angle_to_mirror_normal = angle
|
|---|
| 314 | print "angle:", angle/np.pi*180., "deg"
|
|---|
| 315 |
|
|---|
| 316 | # the rotation matrix for the rotation of *v* around normal_of_reflection_plane is
|
|---|
| 317 | R = make_rotation_matrix( normal_of_reflection_plane, 2*angle )
|
|---|
| 318 |
|
|---|
| 319 | print "R"
|
|---|
| 320 | pprint(R)
|
|---|
| 321 |
|
|---|
| 322 | new_direction = np.dot( R, v)
|
|---|
| 323 | photon.new_direction = new_direction
|
|---|
| 324 |
|
|---|
| 325 | print "old direction", v, length(v)
|
|---|
| 326 | print "new direction", new_direction, length(new_direction)
|
|---|
| 327 | print "mirror center", mirror.pos
|
|---|
| 328 | print "interception point", i
|
|---|
| 329 | print "center of focal plane", focal_plane.pos
|
|---|
| 330 | """
|
|---|
| 331 |
|
|---|
| 332 | # new the photon has a new direction *new_direction* and is starting
|
|---|
| 333 | # from the intersection point *i*
|
|---|
| 334 | # now I want to find out where there focal plane is hit.
|
|---|
| 335 | # So I have to repeat the stuff from up there
|
|---|
| 336 |
|
|---|
| 337 | #print "np.dot(focal_plane.dir,new_direction))", np.dot(focal_plane.dir,new_direction)
|
|---|
| 338 |
|
|---|
| 339 | lambda_1 = (np.dot(focal_plane.dir ,(focal_plane.pos - i)) / np.dot(focal_plane.dir,new_direction))
|
|---|
| 340 |
|
|---|
| 341 | #print "lambda_1", lambda_1
|
|---|
| 342 | focal_plane_spot = lambda_1 * new_direction + i
|
|---|
| 343 | #print "focal_plane_spot",focal_plane_spot
|
|---|
| 344 | photon.hit = True
|
|---|
| 345 | focal_plane_pos = focal_plane_spot - focal_plane.pos
|
|---|
| 346 | photon.focal_plane_pos =focal_plane_pos
|
|---|
| 347 | #photon.hit = True
|
|---|
| 348 | if sh.length(focal_plane_pos) <= focal_plane.size:
|
|---|
| 349 | photon.hit = True
|
|---|
| 350 | else:
|
|---|
| 351 | photon.hit = False
|
|---|
| 352 | return photon
|
|---|
| 353 | # now as a final step we have to find the coordinates of the vector
|
|---|
| 354 | # from the center of the focal plane to the spot where the photon
|
|---|
| 355 | # actually hit the focal plane, as if the plane was not turned.
|
|---|
| 356 | # so if we turn the plane back into the x-y-plane
|
|---|
| 357 | # our *focal_plane_pos* vector has only two coordinates x,y, which are non-zero.
|
|---|
| 358 | # so lets do that.
|
|---|
| 359 | # in order to do so, we need the angles, by which the telescope
|
|---|
| 360 | # was turned .. we hae made them global variables !!ugly i know!!
|
|---|
| 361 |
|
|---|
| 362 | R = make_rotation_matrix( np.array([0,1,0]), -1.*telescope_theta/ 180. *np.pi )
|
|---|
| 363 | turned_focal_plane_pos = np.dot( R, focal_plane_pos)
|
|---|
| 364 | R = make_rotation_matrix( np.array([-1,0,0]), -1.*telescope_phi/ 180. *np.pi )
|
|---|
| 365 | turned_focal_plane_pos = np.dot( R, turned_focal_plane_pos)
|
|---|
| 366 | photon.turned_focal_plane_pos = turned_focal_plane_pos
|
|---|
| 367 | #if np.abs(turned_focal_plane_pos[2] ) > 1e-12:
|
|---|
| 368 | # print turned_focal_plane_pos[2]
|
|---|
| 369 | # raise Exception("the z-coordinate should be zero but is larger than 1e-12")
|
|---|
| 370 |
|
|---|
| 371 |
|
|---|
| 372 |
|
|---|
| 373 | #print "distance from focal plane center=", length(focal_plane_spot-focal_plane.pos)
|
|---|
| 374 | else:
|
|---|
| 375 | photon.hit = False
|
|---|
| 376 | return photon
|
|---|
| 377 |
|
|---|
| 378 |
|
|---|
| 379 | if __name__ == '__main__':
|
|---|
| 380 |
|
|---|
| 381 | # these three things define my telescope today:
|
|---|
| 382 | # * a set off mirrors, read from a file
|
|---|
| 383 | # * a focal plane, which has a postion, a direction and a size
|
|---|
| 384 | # * and a mirror alignmen point, which is needed to construct the mirror
|
|---|
| 385 | # normal vectors.
|
|---|
| 386 | #
|
|---|
| 387 | mirrors = read_reflector_definition_file( "030/fact-reflector.txt" )
|
|---|
| 388 | focal_plane = Focal_Plane( pos=np.array([0.,0.,978.132/2.]), # center of focal_plane
|
|---|
| 389 | dir=np.array([0., 0., 1.]), # direction of view
|
|---|
| 390 | size=20 ) # radius in cm
|
|---|
| 391 |
|
|---|
| 392 | mirror_alignmen_point = Point( pos=np.array([0.,0.,978.132]) )
|
|---|
| 393 |
|
|---|
| 394 |
|
|---|
| 395 | # Now we read the corsika file, which will give us a few thousand photons
|
|---|
| 396 | # to work with. But in order to work with these photons,
|
|---|
| 397 | # we need to turn our telescope into the right direction.
|
|---|
| 398 | # In order to find the right direction, we simply use the mean
|
|---|
| 399 | # direction of the photons in the corsika file.
|
|---|
| 400 | #
|
|---|
| 401 | # In addition we change a little bit in the output format of
|
|---|
| 402 | # the cosika files ... we move all the photons so, they hit a 5m circle
|
|---|
| 403 | # I can't explain that all here, please ask me or wait for the docu :-(
|
|---|
| 404 | #
|
|---|
| 405 | print "working on corsika file: ", sys.argv[1]
|
|---|
| 406 | corsika = readcorsika.read_corsika_file(sys.argv[1])
|
|---|
| 407 |
|
|---|
| 408 | # so first we want to loop over all events in corsika
|
|---|
| 409 | # and move the photons of each event.
|
|---|
| 410 | # in addition we want to find the mean direction of *all* photons in corsika.
|
|---|
| 411 | uv_event_means = []
|
|---|
| 412 | for event in corsika.events:
|
|---|
| 413 | # jump over empty events... I wonder why they exist...
|
|---|
| 414 | if event.info_dict['num_photons'] == 0:
|
|---|
| 415 | continue
|
|---|
| 416 | core_loc = np.array(event.info_dict['core_loc'])
|
|---|
| 417 | # subtract the core location of this event from the x and y coordinates
|
|---|
| 418 | event.photons[:,1:3] -= core_loc
|
|---|
| 419 | uv_mean = event.photons[:,3:5].mean(axis=0)
|
|---|
| 420 | uv_event_means.append(uv_mean)
|
|---|
| 421 | uv_event_means = np.array(uv_event_means)
|
|---|
| 422 |
|
|---|
| 423 | u,v = uv_event_means.mean(axis=0).tolist()
|
|---|
| 424 | print "u,v mean =", u,v
|
|---|
| 425 | theta, phi = sh.uv_to_theta_phi(u,v)
|
|---|
| 426 | theta = theta
|
|---|
| 427 | phi = phi /2.
|
|---|
| 428 | print "theta, phi mean =", theta, phi
|
|---|
| 429 |
|
|---|
| 430 | # turn the telescope
|
|---|
| 431 | # ALARM ... the axis has a minus here .. it works ... but I don't know why.
|
|---|
| 432 | # I turn the telescope here. Around the negative x-axis
|
|---|
| 433 | print mirrors[0].dir
|
|---|
| 434 | turning_axis = np.array([-1,0,0])
|
|---|
| 435 | for mirror in mirrors:
|
|---|
| 436 | mirror.turn( turning_axis, phi)
|
|---|
| 437 | focal_plane.turn( turning_axis, phi)
|
|---|
| 438 | mirror_alignmen_point.turn( turning_axis, phi)
|
|---|
| 439 | # ... and around the y-axis..
|
|---|
| 440 | turning_axis = np.array([0,1,0])
|
|---|
| 441 | for mirror in mirrors:
|
|---|
| 442 | mirror.turn( turning_axis, theta)
|
|---|
| 443 | focal_plane.turn( turning_axis, theta)
|
|---|
| 444 | mirror_alignmen_point.turn( turning_axis, theta)
|
|---|
| 445 | print mirrors[0].dir
|
|---|
| 446 | #
|
|---|
| 447 | # the axes here were found out by trial and error... I still have to find out
|
|---|
| 448 | # which axis is which, in which program and so on. This is pretty confusing still
|
|---|
| 449 | # for me, but on the other hand ... left right up down ... where is the difference :-)
|
|---|
| 450 |
|
|---|
| 451 | global telescope_phi
|
|---|
| 452 | global telescope_theta
|
|---|
| 453 | telescope_phi = phi
|
|---|
| 454 | telescope_theta = theta
|
|---|
| 455 |
|
|---|
| 456 |
|
|---|
| 457 | for event_counter, event in enumerate(corsika.events):
|
|---|
| 458 | print event_counter
|
|---|
| 459 | event.photons_who_hit = []
|
|---|
| 460 | if event.photons is None:
|
|---|
| 461 | continue
|
|---|
| 462 | for photon in event.photons:
|
|---|
| 463 | photon = Photon(photon)
|
|---|
| 464 | photon = reflect_photon( photon, mirrors )
|
|---|
| 465 | if photon.hit:
|
|---|
| 466 | event.photons_who_hit.append(photon)
|
|---|
| 467 | if event_counter > 100:
|
|---|
| 468 | break
|
|---|
| 469 |
|
|---|
| 470 | g = TGraph2D()
|
|---|
| 471 | # g2 = TGraph2D()
|
|---|
| 472 | h = TH2F("h","title",196,-22,22,196,-22,22)
|
|---|
| 473 |
|
|---|
| 474 | graph_point_counter = 0
|
|---|
| 475 | for ev in corsika.events:
|
|---|
| 476 | if not hasattr(ev, "photons_who_hit"):
|
|---|
| 477 | continue
|
|---|
| 478 | for ph in ev.photons_who_hit:
|
|---|
| 479 | tfpp = ph.turned_focal_plane_pos
|
|---|
| 480 | h.Fill(tfpp[0], tfpp[1])
|
|---|
| 481 |
|
|---|
| 482 | fpp = ph.mirror_intersection
|
|---|
| 483 | #fpp = ph.focal_plane_pos
|
|---|
| 484 | g.SetPoint(graph_point_counter, fpp[0],fpp[1],fpp[2])
|
|---|
| 485 | graph_point_counter += 1
|
|---|
| 486 |
|
|---|
| 487 | c1 = TCanvas("c1","c1",0,0,500,500)
|
|---|
| 488 | g.SetMarkerStyle(20)
|
|---|
| 489 | g.Draw("pcol")
|
|---|
| 490 | c1.Update()
|
|---|
| 491 |
|
|---|
| 492 | c2 = TCanvas("c2","c2",0,500,500,500)
|
|---|
| 493 | h.Draw("colz")
|
|---|
| 494 | c2.Update()
|
|---|
