source: fact/tools/pyscripts/simulation/reflector.py@ 17301

Last change on this file since 17301 was 14805, checked in by neise, 12 years ago
initial commit
  • Property svn:executable set to *
File size: 18.8 KB
Line 
1#!/usr/bin/python -itt
2
3import struct
4import sys
5import numpy as np
6from pprint import pprint
7import rlcompleter
8import readline
9readline.parse_and_bind('tab: complete')
10from ROOT import *
11import readcorsika
12import matplotlib.pyplot as plt
13
14from Turnable import *
15
16import simulation_helpers as sh
17
18class 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
46class 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
93class 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
102class Point( Turnable ):
103 """
104 """
105 def __init__(self, pos):
106 self.pos = pos
107 self._turnables = ["pos"]
108
109
110
111def 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
148def 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
379if __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()
Note: See TracBrowser for help on using the repository browser.