#!/usr/bin/python -itt import struct import sys from pprint import pprint import rlcompleter import readline readline.parse_and_bind('tab: complete') #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): def __init__(self, header): self.header = header self.footer = None self.events = [] self.header_dict = {} self.makeHeaderDict() def makeHeaderDict(self): """ this is just moving the plain 272-float numbers from the original header int a dict. Not sure, if the key naming is well done :-/ """ if not len(self.header) == 272: raise TypeException('self.header length must be 272, but is: '+str(len(self.header))) h = self.header d = self.header_dict d['run number']=h[0] d['date of begin run']=int(round(h[1])) d['version of program']=h[2] n_obs_levels = int(round(h[3])) if (n_obs_levels < 1) or (n_obs_levels > 10): ValueException('number of observation levels n must be 0 < n < 11, but is: '+str(h[3])) obs_levels = [] for i in range(4,4+n_obs_levels): obs_levels.append(h[i]) d['observation levels']=tuple(obs_levels) del obs_levels, n_obs_levels d['slope of energy spektrum']=h[14] d['energy range']=(h[15],h[16]) d['flag for EGS4 treatment of em. component'] = h[17] d['flag for NKG treatment of em. component'] = h[18] d['kin. energy cutoff for hadrons in GeV'] = h[19] d['kin. energy cutoff for muons in GeV'] = h[20] d['kin. energy cutoff for electrons in GeV'] = h[21] d['energy cutoff for photons in GeV'] = h[22] physical_constants = [] for i in range(23,73): physical_constants.append(h[i]) d['phyiscal constants']=tuple(physical_constants) del physical_constants d['X-displacement of inclined observation plane'] = h[73] d['Y-displacement of inclined observation plane'] = h[74] d['Z-displacement of inclined observation plane'] = h[75] d['theta angle of normal vector of inclined observation plane'] = h[76] d['phi angle of normal vector of inclined observation plane'] = h[77] # now some constants, I don't understand cka = [] for i in range(93,133): cka.append(h[i]) d['CKA'] = tuple(cka) del cka ceta = [] for i in range(133,138): ceta.append(h[i]) d['CETA'] = tuple(ceta) del ceta cstrba = [] for i in range(138,149): cstrba.append(h[i]) d['CSTRBA'] = tuple(cstrba) del cstrba d['scatter range in x direction for Cherenkov'] = h[246] d['scatter range in y direction for Cherenkov'] = h[247] hlay = [] for i in range(248,253): hlay.append(h[i]) d['HLAY'] = tuple(hlay) del hlay aatm = [] for i in range(253,258): aatm.append(h[i]) d['AATM'] = tuple(aatm) del aatm batm = [] for i in range(258,263): batm.append(h[i]) d['BATM'] = tuple(batm) del batm catm = [] for i in range(263,268): catm.append(h[i]) d['CATM'] = tuple(catm) del catm d['NFLAIN'] = h[268] d['NFLDIF'] = h[269] d['NFLPI0 + 100 x NFLPIF'] = h[270] d['NFLCHE + 100 x NFRAGM'] = h[271] def __repr__(self): print " Run Header " print "="*70 pprint(self.header_dict) print "="*70 print " End of: Run Header " print "="*70 return "" class Event(object): def __init__(self, header): self.header = header self.footer = None self.photons = None self.header_dict = {} self.makeHeaderDict() # Info dict should contain just some quick infor about # this event. # * core location # * energy # * primary incident angle in radian # * # * number of photons self.info_dict = {} def makeInfoDict(self): d = self.info_dict d['energy'] = self.header_dict['total energy in GeV'] d['angle_zenith'] = self.header_dict['angle in radian: (zenith, azimuth)'][0] d['angle_azimuth'] = self.header_dict['angle in radian: (zenith, azimuth)'][1] d['core_loc'] = self.header_dict['core location for scattered events in cm: (x,y)'][0] d['event_id'] = self.header_dict['event number'] d['momentum'] = np.array(self.header_dict['momentum in GeV/c in (x, y, -z) direction;']) #print "self.photons is None", self.photons is None if self.photons is None: d['num_photons'] = 0 else: d['num_photons'] = self.photons.shape[0] def printInfoDict(self): print " Info Dict " print "="*70 pprint(self.info_dict) print "="*70 print " End of: Info Dict " print "="*70 return "" def makeHeaderDict(self): if not len(self.header) == 272: raise TypeException('self.header length must be 272, but is: '+str(len(self.header))) h = self.header d = self.header_dict d['event number'] = long( round(h[0]) ) d['particle id (particle code or A x 100 + Z for nuclei)'] = long( round(h[1]) ) d['total energy in GeV'] = h[2] d['starting altitude in g/cm2'] = h[3] d['number of first target if fixed'] = h[4] d['z coordinate (height) of first interaction in cm'] = h[5] d['momentum in GeV/c in (x, y, -z) direction;'] = (h[6],h[7],h[8]) d['angle in radian: (zenith, azimuth)'] = (h[9],h[10]) n_random_number_sequences = int(round(h[11])) if (n_random_number_sequences < 1) or (n_random_number_sequences > 10): ValueException('number of random number sequences n must be 0 < n < 11, but is: '+str(h[11])) rand_number_sequences = [] for i in range(12,12 + 3*n_random_number_sequences,3): seed = long( round(h[i]) ) number_of_calls = long(round(h[i+2]))*1e6 + long(round(h[i+1])) rand_number_sequences.append( (seed, number_of_calls) ) d['random number sequences: (seed, # of offset random calls)']=tuple(rand_number_sequences) # these are already in run_header ... # can be used for sanity check # ------------------------------------------------------- d['run number'] = h[42] d['date of begin run (yymmdd)'] = int( round(h[43]) ) d['version of program'] = h[44] n_obs_levels = int(round(h[45])) if (n_obs_levels < 1) or (n_obs_levels > 10): ValueException('number of observation levels n must be 0 < n < 11, but is: '+str(h[3])) obs_levels = [] for i in range(46,46+n_obs_levels): obs_levels.append(h[i]) d['observation levels']=tuple(obs_levels) del obs_levels, n_obs_levels d['slope of energy spektrum']=h[56] d['energy range']=(h[57],h[58]) d['kin. energy cutoff for hadrons in GeV'] = h[59] d['kin. energy cutoff for muons in GeV'] = h[60] d['kin. energy cutoff for electrons in GeV'] = h[61] d['energy cutoff for photons in GeV'] = h[62] d['NFLAIN'] = h[63] d['NFLDIF'] = h[64] d['NFLPI0'] = h[65] d['NFLPIF'] = h[66] d['NFLCHE'] = h[67] d['NFRAGM'] = h[68] d["Earth's magnetic field in uT: (x,z)"] = ( h[69], h[70] ) d['flag for activating EGS4'] = h[71] d['flag for activating NKG'] = h[72] d['low-energy hadr. model flag (1.=GHEISHA, 2.=UrQMD, 3.=FLUKA)'] = h[73] d['high-energy hadr. model flag (0.=HDPM,1.=VENUS, 2.=SIBYLL,3.=QGSJET, 4.=DPMJET, 5.=NE X US, 6.=EPOS)'] = h[74] d['CERENKOV Flag (is a bitmap --> usersguide)'] = hex(int(round(h[75]))) d['NEUTRINO flag'] = h[76] d['CURVED flag (0=standard, 2=CURVED)'] = h[77] d['computer flag (3=UNIX, 4=Macintosh)'] = h[78] d['theta interval (in degree): (lower, upper edge) '] = ( h[79], h[80] ) d['phi interval (in degree): (lower, upper edge) '] = ( h[81], h[82] ) d['Cherenkov bunch size in the case of Cherenkov calculations'] = h[83] d['number of Cherenkov detectors in (x, y) direction'] = (h[84], h[85]) d['grid spacing of Cherenkov detectors in cm (x, y) direction'] = ( h[86], h[87]) d['length of each Cherenkov detector in cm in (x, y) direction'] = ( h[88], h[89]) d['Cherenkov output directed to particle output file (= 0.) or Cherenkov output file (= 1.)'] = h[90] d['angle (in rad) between array x-direction and magnetic north'] = h[91] d['flag for additional muon information on particle output file'] = h[92] d['step length factor for multiple scattering step length in EGS4'] = h[93] d['Cherenkov bandwidth in nm: (lower, upper) end'] = ( h[94], h[95] ) d['number i of uses of each Cherenkov event'] = h[96] core_location_for_scattered_events_x_location = [] core_location_for_scattered_events_y_location = [] for i in range(97,117): core_location_for_scattered_events_x_location.append(h[i]) for i in range(117,137): core_location_for_scattered_events_y_location.append(h[i]) d['core location for scattered events in cm: (x,y)'] = zip( core_location_for_scattered_events_x_location, core_location_for_scattered_events_y_location) d['SIBYLL interaction flag (0.= no SIBYLL, 1.=vers.1.6; 2.=vers.2.1)'] = h[137] d['SIBYLL cross-section flag (0.= no SIBYLL, 1.=vers.1.6; 2.=vers.2.1)'] = h[138] d['QGSJET interact. flag (0.=no QGSJET, 1.=QGSJETOLD,2.=QGSJET01c, 3.=QGSJET-II)'] = h[139] d['QGSJET X-sect. flag (0.=no QGSJET, 1.=QGSJETOLD,2.=QGSJET01c, 3.=QGSJET-II)'] = h[140] d['DPMJET interaction flag (0.=no DPMJET, 1.=DPMJET)'] = h[141] d['DPMJET cross-section flag (0.=no DPMJET, 1.=DPMJET)'] = h[142] d['VENUS/NE X US/EPOS cross-section flag (0=neither, 1.=VENUSSIG,2./3.=NEXUSSIG, 4.=EPOSSIG)'] = h[143] d['muon multiple scattering flag (1.=Moliere, 0.=Gauss)'] = h[144] d['NKG radial distribution range in cm'] = h[145] d['EFRCTHN energy fraction of thinning level hadronic'] = h[146] d['EFRCTHN x THINRAT energy fraction of thinning level em-particles'] = h[147] d['actual weight limit WMAX for thinning hadronic'] = h[148] d['actual weight limit WMAX x WEITRAT for thinning em-particles'] = h[149] d['max. radius (in cm) for radial thinning'] = h[150] d['viewing cone VIEWCONE (in deg): (inner, outer) angle'] = (h[151], h[152]) d['transition energy high-energy/low-energy model (in GeV)'] = h[153] d['skimming incidence flag (0.=standard, 1.=skimming)'] = h[154] d['altitude (cm) of horizontal shower axis (skimming incidence)'] = h[155] d['starting height (cm)'] = h[156] d['flag indicating that explicite charm generation is switched on'] = h[156] d['flag for hadron origin of electromagnetic subshower on particle tape'] = h[157] d['flag for observation level curvature (CURVOUT) (0.=flat, 1.=curved)'] = h[166] def printHeaderDict(self): print " Event Header " print "="*70 pprint(self.header_dict) print "="*70 print " End of: Event Header " print "="*70 def __repr__(self): self.printInfoDict() return "" ############################################################################## ############################################################################## # # Two helper functions for plotting the contents of a corsika run follow # def plot_x_y(run, x, y, bins): plt.ion() fig = plt.figure() if y == None: # find min and max of x xmin = np.inf xmax = -np.inf for ev in run.events: if ev.photons[:,x].min() < xmin: xmin = ev.photons[:,x].min() if ev.photons[:,x].max() > xmax: xmax = ev.photons[:,x].max() ra = np.array([xmin,xmax]) Hsum = None for ev in run.events: p = ev.photons H,xb = np.histogram(p[:,x], bins=bins, range=ra) if Hsum == None: Hsum = H else: Hsum += H width = 0.7*(xb[1]-xb[0]) center = (xb[:-1]+xb[1:])/2 plt.bar(center, Hsum, align = 'center', width = width) else: # find min and max of x and y xmin = np.inf xmax = -np.inf ymin = np.inf ymax = -np.inf for ev in run.events: #print x,".min=", ev.photons[:,x].min() #print x,".max=", ev.photons[:,x].max() #print y,".min=", ev.photons[:,y].min() #print y,".max=", ev.photons[:,y].max() if ev.photons is None: continue if ev.photons[:,x].min() < xmin: xmin = ev.photons[:,x].min() if ev.photons[:,x].max() > xmax: xmax = ev.photons[:,x].max() if ev.photons[:,y].min() < ymin: ymin = ev.photons[:,y].min() if ev.photons[:,y].max() > ymax: ymax = ev.photons[:,y].max() #xmin, ymin = -1, -1 #xmax, ymax = 1,1 # gather complete histo ra = np.array([xmin,xmax,ymin,ymax]).reshape(2,2) Hsum = None for ev in run.events: p = ev.photons if p is None: continue H,xb,yb = np.histogram2d(p[:,x], p[:,y],bins=[bins,bins], range=ra) if Hsum == None: Hsum = H else: Hsum += H extent = [ yb[0] , yb[-1] , xb[0] , xb[-1] ] plt.imshow(Hsum, extent=extent, interpolation='nearest') plt.colorbar() def plot_xy_rel_core_location(run, bins): plt.ion() fig = plt.figure() # find min and max of x and y xmin = np.inf xmax = -np.inf ymin = np.inf ymax = -np.inf for ev in run.events: #print x,".min=", ev.photons[:,x].min() #print x,".max=", ev.photons[:,x].max() #print y,".min=", ev.photons[:,y].min() #print y,".max=", ev.photons[:,y].max() if ev.photons is None: continue # get core-location from event header dictionary cl = np.array(ev.header_dict['core location for scattered events in cm: (x,y)'][0]) xy = ev.photons[:,1:3] - cl if xy[:,0].min() < xmin: xmin = xy[:,0].min() if xy[:,0].max() > xmax: xmax = xy[:,0].max() if xy[:,1].min() < ymin: ymin = xy[:,1].min() if xy[:,1].max() > ymax: ymax = xy[:,1].max() #xmin, ymin = -1, -1 #xmax, ymax = 1,1 # gather complete histo ra = np.array([xmin,xmax,ymin,ymax]).reshape(2,2) Hsum = None for index, ev in enumerate(run.events): if ev.photons is None: continue # get core-location from event header dictionary cl = np.array(ev.header_dict['core location for scattered events in cm: (x,y)'][0]) xy = ev.photons[:,1:3] - cl H,xb,yb = np.histogram2d(xy[:,0], xy[:,1], bins=[bins,bins], range=ra) if Hsum == None: Hsum = H else: Hsum += H #print "index: ", index, "H.sum:", H.sum() #print "index: ", index, "Hsum.sum:", Hsum.sum() extent = [ yb[0] , yb[-1] , xb[0] , xb[-1] ] plt.imshow(Hsum, extent=extent, interpolation='nearest') plt.colorbar() return Hsum def read_corsika_file( filename ): f = open(filename, 'rb') while True: # the only break out of this while is finding a run footer block - "RUNE" # Each "block" in a corika file is always 273 "words" long. # a work is a 4-byte value, usually a double, but might also be an int. # The first "word" of a block is sometimes a sequence of 4 ascii # charaters, in order to tell the reader, what type of block follows next. # these sequences can be: # * "RUNH" - a run header block # * "RUNE" - a run footer block # * "EVTH" - an event header block # * "EVTE" - an event footer block # the corsika manual states more of there, but those are not # expected in a corsika "cer-file". block = f.read(273*4) if not block: raise Exception('End of file '+ filename +' reached, before ' + '"RUNE" word found. File might be corrupted.') if not len(block) == 273*4: raise Exception("""While reading %s a block smaller than 273*4 bytes was read, this might be due to non-blocking read which is clearly a bug in this software. you might want to simply try again. """ % (filename,)) # what ever block we have, we treat the first 4 bytes as # a 4-character string. In case it is not one of the # "block descriptors", we assume it is a data-block. block_type = block[0:4] # this reader has no error handling of the underlying file at all. # so in case there is an EVTH found before a RUNH the reader will # probably totally freak out :-) # This is a task for captain state-machine :-) if block_type == 'RUNH': # if we hava a run-header block, the trailing 272 words are # unpacked into a list of floats. from this list a Run-object # is created, which is mainly empty up to that point. run_header = struct.unpack( '272f', block[4:] ) run = Run(run_header) elif block_type == 'EVTH': # directly after the RUNH an event header is expected. # so if that event header is found. again the trailing # 272 words areunpacked into a list of 272 floats and used to # call the constructor of an Event instance. # this event, still it is mainly empty, is appended to the # list of events, maintained by the run object. # # go on reading in the else-block of this of condition. event_header = struct.unpack( '272f', block[4:] ) event = Event(event_header) run.events.append(event) elif block_type == 'EVTE': # after all data blocks for a certain event have been found # we will eventually find a event footer block. # the footer does not contain much information, but it is still # safed with the current event. # It could be used to check, if we found a footer, which fits to the # event header, but that is not done here. event.footer = struct.unpack( '272f', block[4:] ) # When footer is found, we know that the event is # fully read. So we can make a python dict, which contains # some short info about his event, like the number of photons # which is not stored in the event header. event.makeInfoDict() event = None elif block_type == 'RUNE': # when the run footer block is found, we know we have # reached the end of the file, since as far as I have seen # there is only one run in one corsika file. # In addition one can find a lot of zeroes at the end of a corsika # file after the run footer block, so it is not safe to # simply read until the end of file and try to parse everything. run.footer = struct.unpack( '272f', block[4:] ) break else: # since there was no block-descripter found, we assume this # block is a data block. the 273 words of this block are # again unpacked into 273 floats this time. # from the corsika manual we know, that a data block is grouped # into 39 groups of 7 words, so we form an np.array of shape # 39 x 7 data_block = struct.unpack( '273f', block ) data_array = np.array( data_block, dtype=np.float32 ).reshape(39,7) # the data we read in here, does not comply with the original # data sheet. # corsika was edited by somebody calles Dorota from Lodz in order to write the # wavelength of the simulated cherenkov photons and some other info. # I did not find any docu, but had a look into the code, which # luckily is stored right next to the binary. :-) # according to that code the 1st element of the 7-element-tuple # belonging to each cherenkov-photon stores 3 types of information: # a parameter called J ... which might be the parent particles type # a parameter called IMOV, which seems to be constant = 1 # and which seems to have something to do with reusing showers # a parameter called WAVELENGTH which seems to be the cherekov # photon wavelength in nm. # # the code looks like this: # J*100000. + IMOV*1000. + WAVELENGTH new_extra_data = np.zeros( (39,3), dtype=np.float32 ) for index, row in enumerate(data_array): integer = long(round(row[0])) J = integer/100000 IMOV = (integer-J*100000)/1000 WAVELENGTH = np.mod(row[0], 1000.) new_extra_data[index,0] = J new_extra_data[index,1] = IMOV new_extra_data[index,2] = WAVELENGTH data_array = np.concatenate( (data_array,new_extra_data), axis=1) # so at this point we have made from the 39x7 shaped array an # 39x10 shaped array. # the contents of this array are 4 byte floats and have # the following meansings: # d = data_array[0] # d[0] - encoded info from dorota. see index:7,8,9 # d[1] - x position in cm; x-axis points north in corsika # d[2] - y position in cm; y-axis points west in corsika # both positions are measured at the "observation level" aka ground # d[3] - direction cosine to x-axis, called "u" in corsika # d[4] - direction cosine to y-axis, called "v" in corsika # d[5] - time since first interaction in ns # d[6] - height of production in cm; # # now we can't be sure, since we have it reverse engineered. # d[7] - "J", maybe particle ID of ch. light emitting particle # d[8] - "IMOV", apparently always 1, connection with reusing of showers # d[9] - "WAVELENGTH" wavelength of the photon in nm. # since a data block has a fixed size, it can contain # empty rows, we want to find the first empty row here. first_empty_row_index = False for row_index, row in enumerate(data_array): if (row == 0.).all(): first_empty_row_index = row_index break # and here we simply cut off all empty rows. if first_empty_row_index: data_array = data_array[ 0:first_empty_row_index, : ] # here we append our data_array to the np.array, which is stored # in the current event. keep in mind that the current event can simply # be retrieved by "event". if event.photons == None: event.photons = data_array else: event.photons = np.concatenate( (event.photons, data_array), axis=0 ) f.close() run = add_all_photon_table(run) return run if __name__ == '__main__': run = read_corsika_file( sys.argv[1] ) print 'Run object "run" created.' print 'type: print run' print 'or' print 'for event in run.events:' print ' for photon in event.photons:' print ' print photon' print print 'to loop over everyting...' #plot_x_y(run, 3,4,1000) #plt.title("photon distribution ... in theta vs phi ... or other way round") #plot_x_y(run, 1,2,1000) #plt.title("photon distribution ... x vs y .. or other way round") #Hsum = plot_xy_rel_core_location(run, 100) #plt.title("photon distribution 100 bins; file:" + sys.argv[1]) #plt.xlabel('photon y location relative to "core location" in cm') #plt.ylabel('photon x location relative to "core location" in cm')