Changeset 14821 for fact/tools/pyscripts/simulation/readcorsika.py

Timestamp:

02/01/13 08:11:45 (12 years ago)

Author:

neise

Message:

new full photon table added to Run object after reading the entire file.

File:

• fact/tools/pyscripts/simulation/readcorsika.py (modified) (4 diffs)

fact/tools/pyscripts/simulation/readcorsika.py

@@ -3 +3 @@
 import struct
 import sys
-import numpy as np
+
 from pprint import pprint
 
@@ -10 +10 @@
 readline.parse_and_bind('tab: complete')
 
-import matplotlib.pyplot as plt
+#import matplotlib.pyplot as plt
+
+import numpy as np
+
+
+
+
+def add_all_photon_table( run ):
+    ##########################################################################
+    #
+    # At this point I have (hopefully) successfully read the entire file.
+    # But the corsika coordinate system is not quite what I like.
+    # In corsika the x-axis points north
+    # and the y-axis points north.
+    # And so are the direction cosines.
+    # In addition I later found, that I would like to loop over all photons
+    # in a corsika file, without respect to the event in order to delete all 
+    # those which are absorbed or which definitely do not hit the detector,
+    # because they are to far away.
+    # This should reduce computation time a bit...
+    # 
+    # So what I do now, Is setting up a np.ndarray with two dimensions
+    # ( photon_id , parameters ) so I can use nice numpy array operations.
+    #
+    # the former array contained:
+    #        the *photon_definition_array* pda contains:
+    #        pda[0] - encoded info
+    #        pda[1:3] - x,y position in cm. x:north; y:west; z:down
+    #        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]
+    #
+    # 
+    # the new parameters will be:
+    # 0 - x-position in [cm] - x-axis points north
+    # 1 - y-position in [cm] - y-axis points east
+    # 2 - z-position in [cm] - z-axis points up
+    #
+    # 3 - now the direction cosine to x-axis: was formerly u
+    # 4 - now the direction cosine to y-axis: was formerly -v
+    # 5 - now the direction cosine to z-axis: is 1-[3]-[4]
+    #
+    # the next two are in radians:
+    # 6 - theta (angle between photon path and z-axis): computed from 3, 4, and 5
+    # 7 - phi   (angle between photon path projected into x-y-plane and x-axis)
+    #
+    # 8 - height of production in cm
+    # 9 - wavelength in nm
+    # 10- time since 1st interaction in ns
+    # 11- core_location distance in cm
+    # 12- event id
+    # 13- photon id
+    #
+    # in order to find the length of the array I loop ever the events and count the photons
+    #
+    # we directly take care for the so called core_location
+    
+    num_photons = 0
+    for event in run.events:
+        num_photons += event.info_dict['num_photons']
+        if num_photons == 0:
+            continue
+    run.table = np.zeros((num_photons,14))
+    
+    n = 0
+    for event_id, event in enumerate(run.events):
+        this_event_photons = event.info_dict['num_photons']
+        if this_event_photons == 0:
+            continue
+        core_loc = event.info_dict['core_loc']
+        if core_loc is None:
+            print "core_loc is None!:", event_id
+        if event.photons is None:
+            print "event.photons is None", event_id 
+            print "event.photons is None", num_photons
+            
+        shape = (this_event_photons,)
+        # new x(north)   <--- old x(north)
+        run.table[n:n+this_event_photons,0] = event.photons[:,1]-core_loc[0]
+        # new y(east)  <--- old -y(west)
+        run.table[n:n+this_event_photons,1] = -1. * (event.photons[:,2]-core_loc[1])
+        # new z is zero anyway.
+        run.table[n:n+this_event_photons,2] = np.zeros(shape)
+        
+        # new direct. cos to north-axis    <--- old u
+        run.table[n:n+this_event_photons,3] = event.photons[:,3]
+        # new direct. cos to east-axis   <--- old -v
+        run.table[n:n+this_event_photons,4] = -1. * event.photons[:,4]
+
+        # new 3rd direction cosine         <--- old sqrt(1-u^2-v^2)
+        u = run.table[n:n+this_event_photons,3]
+        v = run.table[n:n+this_event_photons,4]
+        run.table[n:n+this_event_photons,5] = np.sqrt(1-u**2-v**2)
+        
+        # new theta
+        run.table[n:n+this_event_photons,6] = np.arccos(run.table[n:n+this_event_photons,5])
+        # new phi
+        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])
+        
+        # production height in cm
+        run.table[n:n+this_event_photons,8] = event.photons[:,6]
+        # wavelength
+        run.table[n:n+this_event_photons,9] = event.photons[:,9]
+        # time since 1st interaction
+        run.table[n:n+this_event_photons,10] = event.photons[:,5]
+        # core_loc distance 
+        run.table[n:n+this_event_photons,11] = np.ones(shape) * np.linalg.norm(core_loc)
+        
+        # event_id
+        run.table[n:n+this_event_photons,12] = np.ones(shape) * event_id
+        # photon_id
+        run.table[n:n+this_event_photons,13] = np.arange(this_event_photons)
+        
+        n += this_event_photons
+    
+    return run
+
+
+
+
 
 class Run(object):
@@ -444 +566 @@
     plt.colorbar()
     return Hsum
+
+
+
 
 
@@ -609 +734 @@
     
     f.close()
+        
+    run = add_all_photon_table(run)
 
     return run