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