#!/usr/bin/python -itt import struct import sys import numpy as np from pprint import pprint import rlcompleter import readline readline.parse_and_bind('tab: complete') from ROOT import * import readcorsika import matplotlib.pyplot as plt def dot( v1, v2): return (v1*v2).sum() def length( v ): return np.sqrt((v**2).sum()) def cross(v1 , v2, normalize=True): vout = np.zeros( v1.shape ) vout[0] = v1[1]*v2[2] - v1[2]*v2[1] vout[1] = -(v1[0]*v2[2] - v1[2]*v2[0]) vout[2] = v1[0]*v2[1] - v1[1]*v2[0] if normalize: vout /= length(vout) return vout def matrix_times_vector( m , v): vout = v.copy() for index,line in enumerate(m): vout[index] = dot(line,v) return vout def make_rotation_matrix( nn, a ): R = np.zeros( (3,3) ) R[0,0] = nn[0]*nn[0] * (1-np.cos(a)) + np.cos(a) R[1,0] = nn[0]*nn[1] * (1-np.cos(a)) + nn[2]*np.sin(a) R[2,0] = nn[0]*nn[2] * (1-np.cos(a)) - nn[1]*np.sin(a) R[0,1] = nn[0]*nn[1] * (1-np.cos(a)) - nn[2]*np.sin(a) R[1,1] = nn[1]*nn[1] * (1-np.cos(a)) + np.cos(a) R[2,1] = nn[1]*nn[2] * (1-np.cos(a)) + nn[0]*np.sin(a) R[0,2] = nn[0]*nn[2] * (1-np.cos(a)) + nn[1]*np.sin(a) R[1,2] = nn[1]*nn[2] * (1-np.cos(a)) - nn[0]*np.sin(a) R[2,2] = nn[2]*nn[2] * (1-np.cos(a)) + np.cos(a) return R class Thing( object ): """ Thing is just a container for the postion and the direction of something. A Thing can be a particle, or photon or something like that Or it can be a plane, like a mirror-plane or a focal-plane. Then the *dir* vector is the normal vector of the plane, and *pos* is one (possibly important) point inside the plane. """ def __init__(self, pos, dir): self.pos = pos self.dir = dir def turn(self, axis, angle): """ axis might not be normalized and angle might be in degree """ if length(axis) != 1.: axis /= length(axis) angle = angle / 180. *np.pi R = make_rotation_matrix( turning_axis, angle ) self.pos = matrix_times_vector( R, self.pos) self.dir = matrix_times_vector( R, self.dir) def _old_turn( self, theta, phi): theta = theta/180.*np.pi phi = phi/180.*np.pi x= self.pos[0] y= self.pos[1] z= self.pos[2] vx= self.dir[0] vy= self.dir[1] vz= self.dir[2] #print vx,vy,vz #transform into spehrical coordinates print x,y,z r = length(self.pos) p = np.arctan2(y,x) t = np.arccos(z/r) print r,t,p v = length(self.dir) vp = np.arctan2(vy,vx) vt = np.arccos(vz/v) #print v,vt,vp # actual turning takes place t += theta p += phi vt += theta vp += phi #print v,vt,vp #back transform x = r * np.sin(t) * np.cos(p) y = r * np.sin(t) * np.sin(p) z = r * np.cos(t) vx = v * np.sin(vt) * np.cos(vp) vy = v * np.sin(vt) * np.sin(vp) vz = v * np.cos(vt) #print vx,vy,vz # set internal vars self.pos = np.array([x,y,z]) self.dir = np.array([vx,vy,vz]) def __repr__( self ): return "%s(%r)" % (self.__class__, self.__dict__) class Mirror( Thing ): def __init__(self, index, pos, normal_vector, focal_length, hex_size ): super(Mirror,self).__init__(pos, normal_vector) self.index = index self.focal_length = focal_length self.hex_size = hex_size def __repr__( self ): return "%s(%r)" % (self.__class__, self.__dict__) def read_reflector_definition_file( filename ): """ """ mirrors = [] f = open( filename ) for index, line in enumerate(f): if line[0] == '#': continue line = line.split() if len(line) < 8: continue #print line # first 3 colums in the file are x,y,z coordinates of the center # of this mirror in cm, I guess pos = np.array(map(float, line[0:3])) # the next 3 elements are the elements of the normal vector # should be normalized already, so the unit is of no importance. normal_vector = np.array(map(float,line[3:6])) # focal length of this mirror in mm focal_length = float(line[6]) # size of the hexagonal shaped facette mirror, measured as the radius # of the hexagons *inner* circle. hex_size = float(line[8]) mirror = Mirror( index, pos, normal_vector, focal_length, hex_size ) mirrors.append(mirror) return mirrors def reflect_photon( photon, mirrors): """ finds out: which mirror is hit by photon and where and in which angle relative to mirror """ # the line defined by the photon is used to find the intersection point # with the plane of each facette mirror. Then I check, # if the intersection point lies within the limits of the facette mirrors # hexagonal boundaries. # If this is the case I have found the mirror, which is hit, and # can calculate: # the distance of the intersection point from the center of the facette mirror # and the angle relative to the mirror (normal or plane not sure yet) for mirror in mirrors: #facette mirror plane, defined as n . x = d1 . n n = mirror.dir d1 = mirror.pos # line of photon defined as r = lambda * v + d2 v = photon.dir d2 = photon.pos # the intersection coordinates are found by solving # n . (lambda * v + d2) - d1 . n == 0, for lambda=lambda_0 # and then the intersection is: i = lambda_0 * v + d2 # # putting int in another form: # solve: # lambda * n.v + n.d2 - n.d1 == 0 # or # lambda_0 = n.(d1-d2) / n.v # FIXME: if one of the two dot-products is very small, # we shuold take special care maybe # if n.(d1-d2) is very small, this means that the starting point of # the photon is already nearly in the plane, so lambda_0 is expected to # be very small ... erm .. maybe this is actually not a special case # but very good. # of n.v is very small, this means the patch of the photon is nearly # parallel to the plane, so the error ob lambda_0 might be very large. # in addition, this might just tell us, that the mirror is hit under # strange circumstances ... so its not a good candidate and we can already go on. lambda_0 = (dot(n,(d1-d2)) / dot(n,v)) #intersection between line and plane i = lambda_0 * v + d2 # I want the distance beween i and d1 so I can already from the distance find # out if this is our candidate. distance = np.sqrt(((i-d1)*(i-d1)).sum()) #print "photon pos:", d2, "\t dir:", v/length(v) #print "mirror pos:", d1, "\t dir:", n/length(n) #print "lambda_0", lambda_0 #print "intersection :", i #print "distance:",distance if distance <= mirror.hex_size/2.: break else: mirror = None if not mirror is None: photon.mirror_index = mirror.index photon.mirror_intersection = i photon.mirror_center_distance = distance #print "mirror found:", mirror.index , #print "distance", distance # now I have to find out, if the photon is not only in the # right distance but actually has hit the mirror. # this I do like this # i-d1 is a vector in the mirror plane pointing from d1 to the intersection point i. # if I know turn the entire mirror plane so it lies withing the x-y-plane # by applying a simple turning-matrix, then each vector inside the plane with turn into # a nice x,x vector. # now I assume, that the hexagon is "pointing" lets say to into y direction # so I can e.g. say: # x has to be between -30.3 and +30.3 and y has to be # between 35 - m * |x| and -35 + m * |x| ... pretty simple. # maybe one can leave the turning aside, but I like that I can imagine it very nicely # # # I don't do this yet .. since I don't know by heart how a turning matrix looks :-) # so I just simulate round mirrors ###################################################################### # next step, since I know the intersection point, is the new direction. # So I need the normal of the mirror in the intersection point. # since the normal of every mirror is alway pointing to the camera center # this is not difficult. normal_at_intersection = (mirror_alignmen_point.pos - i) / length(mirror_alignmen_point.pos - i) #print "normal_at_intersection",normal_at_intersection angle = np.arccos(dot( v, normal_at_intersection) / (length(v) * length(normal_at_intersection))) photon.angle_to_mirror_normal = angle #print "angle:", angle/np.pi*180., "deg" # okay, now I have the intersection *i*, # the old direction of the photon *v* # and the normalvector at the intersection. ###################################################################### ###################################################################### # I do this now differently. # I will mirror the "point" at the tip of *v* at the line created by # the normalvector at the intersection and the intersection. # this will gibe me a mirrored_point *mp* and the vector from *i* to *mp* # is the *new_direction* it should even be normalized. # 1. step: create plane through the "tip" of *v* and the normal_at_intersection. # 2. step: find crossingpoint on line through *i* and the normal_at_intersection, # 3. step: vector from "tip" of *v* to crossingpoint times 2 points to # the "tip" of *mirrored_v* # plane: n_plane_3 . r = p_plane_3 . n_plane_3 # p_plane_3 = i+v # n_plane_3 = normal_at_intersection # line: r = lambda_3 * v_line_3 + p_line_3 # p_line_3 = i # v_line_3 = normal_at_intersection # create crossing: n_plane_3 . (lambda_3 * v_line_3 + p_line_3) = p_plane_3 . n_plane_3 # <=> lambda_3 = (p_plane_3 - p_line_3 ).n_plane_3 / n_plane_3 . v_line_3 # <=> lambda_3 = (i+v - i).normal_at_intersection / normal_at_intersection . normal_at_intersection # <=> lambda_3 = v.normal_at_intersection lambda_3 = dot(v, normal_at_intersection) #print "lambda_3", lambda_3 crossing_point_3 = lambda_3 * normal_at_intersection + i #print "crossing_point_3", crossing_point_3 from_tip_of_v_to_crossing_point_3 = crossing_point_3 - (i+v) tip_of_mirrored_v = i+v+ 2*from_tip_of_v_to_crossing_point_3 new_direction = tip_of_mirrored_v - i #print "new_direction",new_direction #print "old direction", v photon.new_direction = new_direction ###################################################################### ###################################################################### """ # both directions form a plane, and when I turn the old *v* by # twice the angle between *v* and *normal_at_intersection* # inside this plane then I get the new direction of the photon. # so lets first get the normal of the reflection plane normal_of_reflection_plane =cross( v, normal_at_intersection) print length(normal_of_reflection_plane), "should be one" print length(v), "should be one" print length(normal_at_intersection), "should be one" print dot(v, normal_at_intersection), "should *NOT* be zero" print dot(v, normal_of_reflection_plane), "should be zero" print dot(normal_at_intersection, normal_of_reflection_plane), "should be zero" angle = np.arccos(dot( v, normal_at_intersection) / (length(v) * length(normal_at_intersection))) photon.angle_to_mirror_normal = angle print "angle:", angle/np.pi*180., "deg" # the rotation matrix for the rotation of *v* around normal_of_reflection_plane is R = make_rotation_matrix( normal_of_reflection_plane, 2*angle ) print "R" pprint(R) new_direction = matrix_times_vector( R, v) photon.new_direction = new_direction print "old direction", v, length(v) print "new direction", new_direction, length(new_direction) print "mirror center", mirror.pos print "interception point", i print "center of focal plane", focal_plane.pos """ # new the photon has a new direction *new_direction* and is starting # from the intersection point *i* # now I want to find out where there focal plane is hit. # So I have to repeat the stuff from up there #print "dot(focal_plane.dir,new_direction))", dot(focal_plane.dir,new_direction) lambda_1 = (dot(focal_plane.dir ,(focal_plane.pos - i)) / dot(focal_plane.dir,new_direction)) #print "lambda_1", lambda_1 focal_plane_spot = lambda_1 * new_direction + i #print "focal_plane_spot",focal_plane_spot photon.hit = True photon.focal_plane_pos = focal_plane_spot - focal_plane.pos #print "distance from focal plane center=", length(focal_plane_spot-focal_plane.pos) else: photon.hit = False return photon class Photon( Thing ): """ a photon has not only the direction and position, which a Thing has. but it also has a wavelength and a "time" and a "mother_particle_id" the photon constructor understands the 10-element 1D-np.array which is stored inside a run.event.photons 2D-np.array """ def __init__(self, photon_definition_array ): """ the *photon_definition_array* pda contains: pda[0] - encoded info pda[1:3] - x,y position in cm pda[3:5] - u,v cosines to x,y axis --> so called direction cosines pda[5] - time since first interaction [ns] pda[6] - height of production in cm pda[7] - j ?? pda[8] - imov ?? pda[9] - wavelength [nm] """ pda = photon_definition_array pos = np.array([pda[1],pda[2],0.]) dir = np.array([pda[3],pda[4], np.sqrt(1-pda[3]**2-pda[4]**2) ]) super(Photon,self).__init__(pos, dir) self.wavelength = pda[9] self.time = pda[5] def __repr__( self ): return "%s(%r)" % (self.__class__, self.__dict__) if __name__ == '__main__': mirrors = read_reflector_definition_file( "030/reflector_test_ray.py" ) focal_plane = Thing( pos=np.array([0.,0.,978.132/2.]), # center of focal_plane dir=np.array([0., 0., 1.]) ) # direction of view focal_plane.size = 20 # diameter in cm mirror_alignmen_point = Thing( pos=np.array([0.,0.,978.132]), # center of focal_plane dir=np.array([0., 0., 1.]) ) # direction of view # turn the telescope turning_axis = np.array([-1,0,0]) angle = 30. for mirror in mirrors: mirror.turn( turning_axis, angle) focal_plane.turn( turning_axis, angle) mirror_alignmen_point.turn( turning_axis, angle) #run = readcorsika.read_corsika_file("cer") li=[] photons_who_hit = [] for photon_index in range(2000): if photon_index % 10 == 0: print photon_index # make test photon directly from up to down p = 87.5 + (photon_index/2000.*5.) p = p/180.*np.pi t = 27.5 + (photon_index/2000.*5.) t = t/180.*np.pi dir = np.array([np.sin(t)*np.cos(p),np.sin(t)*np.sin(p),np.cos(t)]) pos = (np.random.rand(3)-0.5)*450 pos[2] = 0. photon = Thing( pos=pos, dir=dir) photon = reflect_photon( photon, mirrors ) if photon.hit: li.append( photon.angle_to_mirror_normal / np.pi * 180.) photons_who_hit.append(photon) #plt.ion() #fig1 = plt.figure() #fig2 = plt.figure() #plt.hist( np.array(li) , bins=100) plt.hold(False) g = TGraph2D() h = TH2F("h","title",100,-100.5,99.5,100,-100.5,99.5) c1 = TCanvas("c1","c1 title",0,0,400,400) c2 = TCanvas("c2","c2 title",0,400,400,400) c = 0 ground_dirs = [] focal_positions = [] print len(photons_who_hit) for ph in photons_who_hit: mi = ph.mirror_intersection #new = Thing( pos=ph.focal_plane_pos, dir = ph.new_direction) #new.turn( np.array([0,1,0]), 0 ) fpp = ph.focal_plane_pos #print mi h.Fill(fpp[0], fpp[1]) g.SetPoint( c, mi[0], mi[1], mi[2]) c += 1 #ground_dirs.append(ph.dir) #focal_positions.append(new.pos) #ground_dirs = np.array(ground_dirs) #focal_positions = np.array(focal_positions) #print ground_dirs.shape, focal_positions.shape #plt.figure(fig1.number) #plt.plot( ground_dirs[:,0], ground_dirs[:,1], '.') #plt.title("ground directions") #plt.figure(fig2.number) #plt.plot( focal_positions[:,0], focal_positions[:,1], '.') #plt.title("focal positions") #raw_input() g.SetMarkerStyle(20) c1.cd() g.Draw("pcol") c1.Update() c2.cd() h.Draw("colz") c2.Update() """ for m in mirrors: mi = m.pos #mi = ph.focal_plane_pos g.SetPoint( c, mi[0], mi[1], mi[2]) c += 1 """