Changeset 14821
- Timestamp:
- 02/01/13 08:11:45 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
fact/tools/pyscripts/simulation/readcorsika.py
r14805 r14821 3 3 import struct 4 4 import sys 5 import numpy as np 5 6 6 from pprint import pprint 7 7 … … 10 10 readline.parse_and_bind('tab: complete') 11 11 12 import matplotlib.pyplot as plt 12 #import matplotlib.pyplot as plt 13 14 import numpy as np 15 16 17 18 19 def add_all_photon_table( run ): 20 ########################################################################## 21 # 22 # At this point I have (hopefully) successfully read the entire file. 23 # But the corsika coordinate system is not quite what I like. 24 # In corsika the x-axis points north 25 # and the y-axis points north. 26 # And so are the direction cosines. 27 # In addition I later found, that I would like to loop over all photons 28 # in a corsika file, without respect to the event in order to delete all 29 # those which are absorbed or which definitely do not hit the detector, 30 # because they are to far away. 31 # This should reduce computation time a bit... 32 # 33 # So what I do now, Is setting up a np.ndarray with two dimensions 34 # ( photon_id , parameters ) so I can use nice numpy array operations. 35 # 36 # the former array contained: 37 # the *photon_definition_array* pda contains: 38 # pda[0] - encoded info 39 # pda[1:3] - x,y position in cm. x:north; y:west; z:down 40 # pda[3:5] - u,v cosines to x,y axis --> so called direction cosines 41 # pda[5] - time since first interaction [ns] 42 # pda[6] - height of production in cm 43 # pda[7] - j ?? 44 # pda[8] - imov ?? 45 # pda[9] - wavelength [nm] 46 # 47 # 48 # the new parameters will be: 49 # 0 - x-position in [cm] - x-axis points north 50 # 1 - y-position in [cm] - y-axis points east 51 # 2 - z-position in [cm] - z-axis points up 52 # 53 # 3 - now the direction cosine to x-axis: was formerly u 54 # 4 - now the direction cosine to y-axis: was formerly -v 55 # 5 - now the direction cosine to z-axis: is 1-[3]-[4] 56 # 57 # the next two are in radians: 58 # 6 - theta (angle between photon path and z-axis): computed from 3, 4, and 5 59 # 7 - phi (angle between photon path projected into x-y-plane and x-axis) 60 # 61 # 8 - height of production in cm 62 # 9 - wavelength in nm 63 # 10- time since 1st interaction in ns 64 # 11- core_location distance in cm 65 # 12- event id 66 # 13- photon id 67 # 68 # in order to find the length of the array I loop ever the events and count the photons 69 # 70 # we directly take care for the so called core_location 71 72 num_photons = 0 73 for event in run.events: 74 num_photons += event.info_dict['num_photons'] 75 if num_photons == 0: 76 continue 77 run.table = np.zeros((num_photons,14)) 78 79 n = 0 80 for event_id, event in enumerate(run.events): 81 this_event_photons = event.info_dict['num_photons'] 82 if this_event_photons == 0: 83 continue 84 core_loc = event.info_dict['core_loc'] 85 if core_loc is None: 86 print "core_loc is None!:", event_id 87 if event.photons is None: 88 print "event.photons is None", event_id 89 print "event.photons is None", num_photons 90 91 shape = (this_event_photons,) 92 # new x(north) <--- old x(north) 93 run.table[n:n+this_event_photons,0] = event.photons[:,1]-core_loc[0] 94 # new y(east) <--- old -y(west) 95 run.table[n:n+this_event_photons,1] = -1. * (event.photons[:,2]-core_loc[1]) 96 # new z is zero anyway. 97 run.table[n:n+this_event_photons,2] = np.zeros(shape) 98 99 # new direct. cos to north-axis <--- old u 100 run.table[n:n+this_event_photons,3] = event.photons[:,3] 101 # new direct. cos to east-axis <--- old -v 102 run.table[n:n+this_event_photons,4] = -1. * event.photons[:,4] 103 104 # new 3rd direction cosine <--- old sqrt(1-u^2-v^2) 105 u = run.table[n:n+this_event_photons,3] 106 v = run.table[n:n+this_event_photons,4] 107 run.table[n:n+this_event_photons,5] = np.sqrt(1-u**2-v**2) 108 109 # new theta 110 run.table[n:n+this_event_photons,6] = np.arccos(run.table[n:n+this_event_photons,5]) 111 # new phi 112 run.table[n:n+this_event_photons,7] = np.arctan2( run.table[n:n+this_event_photons,4], run.table[n:n+this_event_photons,3]) 113 114 # production height in cm 115 run.table[n:n+this_event_photons,8] = event.photons[:,6] 116 # wavelength 117 run.table[n:n+this_event_photons,9] = event.photons[:,9] 118 # time since 1st interaction 119 run.table[n:n+this_event_photons,10] = event.photons[:,5] 120 # core_loc distance 121 run.table[n:n+this_event_photons,11] = np.ones(shape) * np.linalg.norm(core_loc) 122 123 # event_id 124 run.table[n:n+this_event_photons,12] = np.ones(shape) * event_id 125 # photon_id 126 run.table[n:n+this_event_photons,13] = np.arange(this_event_photons) 127 128 n += this_event_photons 129 130 return run 131 132 133 134 13 135 14 136 class Run(object): … … 444 566 plt.colorbar() 445 567 return Hsum 568 569 570 446 571 447 572 … … 609 734 610 735 f.close() 736 737 run = add_all_photon_table(run) 611 738 612 739 return run
Note:
See TracChangeset
for help on using the changeset viewer.