| 1 | #!/usr/bin/python -itt
|
|---|
| 2 |
|
|---|
| 3 | import struct
|
|---|
| 4 | import sys
|
|---|
| 5 |
|
|---|
| 6 | from pprint import pprint
|
|---|
| 7 |
|
|---|
| 8 | import rlcompleter
|
|---|
| 9 | import readline
|
|---|
| 10 | readline.parse_and_bind('tab: complete')
|
|---|
| 11 |
|
|---|
| 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 |
|
|---|
| 135 |
|
|---|
| 136 | class Run(object):
|
|---|
| 137 | def __init__(self, header):
|
|---|
| 138 | self.header = header
|
|---|
| 139 | self.footer = None
|
|---|
| 140 | self.events = []
|
|---|
| 141 | self.header_dict = {}
|
|---|
| 142 | self.makeHeaderDict()
|
|---|
| 143 |
|
|---|
| 144 | def makeHeaderDict(self):
|
|---|
| 145 | """ this is just moving the plain 272-float numbers from the original
|
|---|
| 146 | header int a dict.
|
|---|
| 147 | Not sure, if the key naming is well done :-/
|
|---|
| 148 | """
|
|---|
| 149 | if not len(self.header) == 272:
|
|---|
| 150 | raise TypeException('self.header length must be 272, but is: '+str(len(self.header)))
|
|---|
| 151 |
|
|---|
| 152 | h = self.header
|
|---|
| 153 | d = self.header_dict
|
|---|
| 154 |
|
|---|
| 155 | d['run number']=h[0]
|
|---|
| 156 | d['date of begin run']=int(round(h[1]))
|
|---|
| 157 | d['version of program']=h[2]
|
|---|
| 158 |
|
|---|
| 159 | n_obs_levels = int(round(h[3]))
|
|---|
| 160 | if (n_obs_levels < 1) or (n_obs_levels > 10):
|
|---|
| 161 | ValueException('number of observation levels n must be 0 < n < 11, but is: '+str(h[3]))
|
|---|
| 162 | obs_levels = []
|
|---|
| 163 | for i in range(4,4+n_obs_levels):
|
|---|
| 164 | obs_levels.append(h[i])
|
|---|
| 165 | d['observation levels']=tuple(obs_levels)
|
|---|
| 166 | del obs_levels, n_obs_levels
|
|---|
| 167 |
|
|---|
| 168 | d['slope of energy spektrum']=h[14]
|
|---|
| 169 | d['energy range']=(h[15],h[16])
|
|---|
| 170 | d['flag for EGS4 treatment of em. component'] = h[17]
|
|---|
| 171 | d['flag for NKG treatment of em. component'] = h[18]
|
|---|
| 172 |
|
|---|
| 173 | d['kin. energy cutoff for hadrons in GeV'] = h[19]
|
|---|
| 174 | d['kin. energy cutoff for muons in GeV'] = h[20]
|
|---|
| 175 | d['kin. energy cutoff for electrons in GeV'] = h[21]
|
|---|
| 176 | d['energy cutoff for photons in GeV'] = h[22]
|
|---|
| 177 |
|
|---|
| 178 | physical_constants = []
|
|---|
| 179 | for i in range(23,73):
|
|---|
| 180 | physical_constants.append(h[i])
|
|---|
| 181 | d['phyiscal constants']=tuple(physical_constants)
|
|---|
| 182 | del physical_constants
|
|---|
| 183 |
|
|---|
| 184 | d['X-displacement of inclined observation plane'] = h[73]
|
|---|
| 185 | d['Y-displacement of inclined observation plane'] = h[74]
|
|---|
| 186 | d['Z-displacement of inclined observation plane'] = h[75]
|
|---|
| 187 | d['theta angle of normal vector of inclined observation plane'] = h[76]
|
|---|
| 188 | d['phi angle of normal vector of inclined observation plane'] = h[77]
|
|---|
| 189 |
|
|---|
| 190 | # now some constants, I don't understand
|
|---|
| 191 | cka = []
|
|---|
| 192 | for i in range(93,133):
|
|---|
| 193 | cka.append(h[i])
|
|---|
| 194 | d['CKA'] = tuple(cka)
|
|---|
| 195 | del cka
|
|---|
| 196 |
|
|---|
| 197 | ceta = []
|
|---|
| 198 | for i in range(133,138):
|
|---|
| 199 | ceta.append(h[i])
|
|---|
| 200 | d['CETA'] = tuple(ceta)
|
|---|
| 201 | del ceta
|
|---|
| 202 |
|
|---|
| 203 | cstrba = []
|
|---|
| 204 | for i in range(138,149):
|
|---|
| 205 | cstrba.append(h[i])
|
|---|
| 206 | d['CSTRBA'] = tuple(cstrba)
|
|---|
| 207 | del cstrba
|
|---|
| 208 |
|
|---|
| 209 | d['scatter range in x direction for Cherenkov'] = h[246]
|
|---|
| 210 | d['scatter range in y direction for Cherenkov'] = h[247]
|
|---|
| 211 |
|
|---|
| 212 | hlay = []
|
|---|
| 213 | for i in range(248,253):
|
|---|
| 214 | hlay.append(h[i])
|
|---|
| 215 | d['HLAY'] = tuple(hlay)
|
|---|
| 216 | del hlay
|
|---|
| 217 |
|
|---|
| 218 | aatm = []
|
|---|
| 219 | for i in range(253,258):
|
|---|
| 220 | aatm.append(h[i])
|
|---|
| 221 | d['AATM'] = tuple(aatm)
|
|---|
| 222 | del aatm
|
|---|
| 223 |
|
|---|
| 224 | batm = []
|
|---|
| 225 | for i in range(258,263):
|
|---|
| 226 | batm.append(h[i])
|
|---|
| 227 | d['BATM'] = tuple(batm)
|
|---|
| 228 | del batm
|
|---|
| 229 |
|
|---|
| 230 | catm = []
|
|---|
| 231 | for i in range(263,268):
|
|---|
| 232 | catm.append(h[i])
|
|---|
| 233 | d['CATM'] = tuple(catm)
|
|---|
| 234 | del catm
|
|---|
| 235 |
|
|---|
| 236 | d['NFLAIN'] = h[268]
|
|---|
| 237 | d['NFLDIF'] = h[269]
|
|---|
| 238 | d['NFLPI0 + 100 x NFLPIF'] = h[270]
|
|---|
| 239 | d['NFLCHE + 100 x NFRAGM'] = h[271]
|
|---|
| 240 |
|
|---|
| 241 |
|
|---|
| 242 | def __repr__(self):
|
|---|
| 243 | print " Run Header "
|
|---|
| 244 | print "="*70
|
|---|
| 245 | pprint(self.header_dict)
|
|---|
| 246 | print "="*70
|
|---|
| 247 | print " End of: Run Header "
|
|---|
| 248 | print "="*70
|
|---|
| 249 | return ""
|
|---|
| 250 |
|
|---|
| 251 |
|
|---|
| 252 |
|
|---|
| 253 | class Event(object):
|
|---|
| 254 | def __init__(self, header):
|
|---|
| 255 | self.header = header
|
|---|
| 256 | self.footer = None
|
|---|
| 257 | self.photons = None
|
|---|
| 258 | self.header_dict = {}
|
|---|
| 259 | self.makeHeaderDict()
|
|---|
| 260 |
|
|---|
| 261 | # Info dict should contain just some quick infor about
|
|---|
| 262 | # this event.
|
|---|
| 263 | # * core location
|
|---|
| 264 | # * energy
|
|---|
| 265 | # * primary incident angle in radian
|
|---|
| 266 | # *
|
|---|
| 267 | # * number of photons
|
|---|
| 268 | self.info_dict = {}
|
|---|
| 269 |
|
|---|
| 270 |
|
|---|
| 271 | def makeInfoDict(self):
|
|---|
| 272 |
|
|---|
| 273 | d = self.info_dict
|
|---|
| 274 |
|
|---|
| 275 | d['energy'] = self.header_dict['total energy in GeV']
|
|---|
| 276 | d['angle_zenith'] = self.header_dict['angle in radian: (zenith, azimuth)'][0]
|
|---|
| 277 | d['angle_azimuth'] = self.header_dict['angle in radian: (zenith, azimuth)'][1]
|
|---|
| 278 | d['core_loc'] = self.header_dict['core location for scattered events in cm: (x,y)'][0]
|
|---|
| 279 | d['event_id'] = self.header_dict['event number']
|
|---|
| 280 | d['momentum'] = np.array(self.header_dict['momentum in GeV/c in (x, y, -z) direction;'])
|
|---|
| 281 |
|
|---|
| 282 | #print "self.photons is None", self.photons is None
|
|---|
| 283 | if self.photons is None:
|
|---|
| 284 | d['num_photons'] = 0
|
|---|
| 285 | else:
|
|---|
| 286 | d['num_photons'] = self.photons.shape[0]
|
|---|
| 287 |
|
|---|
| 288 | def printInfoDict(self):
|
|---|
| 289 | print " Info Dict "
|
|---|
| 290 | print "="*70
|
|---|
| 291 | pprint(self.info_dict)
|
|---|
| 292 | print "="*70
|
|---|
| 293 | print " End of: Info Dict "
|
|---|
| 294 | print "="*70
|
|---|
| 295 | return ""
|
|---|
| 296 |
|
|---|
| 297 |
|
|---|
| 298 | def makeHeaderDict(self):
|
|---|
| 299 | if not len(self.header) == 272:
|
|---|
| 300 | raise TypeException('self.header length must be 272, but is: '+str(len(self.header)))
|
|---|
| 301 |
|
|---|
| 302 | h = self.header
|
|---|
| 303 | d = self.header_dict
|
|---|
| 304 |
|
|---|
| 305 | d['event number'] = long( round(h[0]) )
|
|---|
| 306 | d['particle id (particle code or A x 100 + Z for nuclei)'] = long( round(h[1]) )
|
|---|
| 307 |
|
|---|
| 308 | d['total energy in GeV'] = h[2]
|
|---|
| 309 | d['starting altitude in g/cm2'] = h[3]
|
|---|
| 310 | d['number of first target if fixed'] = h[4]
|
|---|
| 311 | d['z coordinate (height) of first interaction in cm'] = h[5]
|
|---|
| 312 | d['momentum in GeV/c in (x, y, -z) direction;'] = (h[6],h[7],h[8])
|
|---|
| 313 | d['angle in radian: (zenith, azimuth)'] = (h[9],h[10])
|
|---|
| 314 |
|
|---|
| 315 | n_random_number_sequences = int(round(h[11]))
|
|---|
| 316 | if (n_random_number_sequences < 1) or (n_random_number_sequences > 10):
|
|---|
| 317 | ValueException('number of random number sequences n must be 0 < n < 11, but is: '+str(h[11]))
|
|---|
| 318 | rand_number_sequences = []
|
|---|
| 319 | for i in range(12,12 + 3*n_random_number_sequences,3):
|
|---|
| 320 | seed = long( round(h[i]) )
|
|---|
| 321 | number_of_calls = long(round(h[i+2]))*1e6 + long(round(h[i+1]))
|
|---|
| 322 | rand_number_sequences.append( (seed, number_of_calls) )
|
|---|
| 323 | d['random number sequences: (seed, # of offset random calls)']=tuple(rand_number_sequences)
|
|---|
| 324 |
|
|---|
| 325 | # these are already in run_header ...
|
|---|
| 326 | # can be used for sanity check
|
|---|
| 327 | # -------------------------------------------------------
|
|---|
| 328 | d['run number'] = h[42]
|
|---|
| 329 | d['date of begin run (yymmdd)'] = int( round(h[43]) )
|
|---|
| 330 | d['version of program'] = h[44]
|
|---|
| 331 |
|
|---|
| 332 | n_obs_levels = int(round(h[45]))
|
|---|
| 333 | if (n_obs_levels < 1) or (n_obs_levels > 10):
|
|---|
| 334 | ValueException('number of observation levels n must be 0 < n < 11, but is: '+str(h[3]))
|
|---|
| 335 | obs_levels = []
|
|---|
| 336 | for i in range(46,46+n_obs_levels):
|
|---|
| 337 | obs_levels.append(h[i])
|
|---|
| 338 | d['observation levels']=tuple(obs_levels)
|
|---|
| 339 | del obs_levels, n_obs_levels
|
|---|
| 340 |
|
|---|
| 341 | d['slope of energy spektrum']=h[56]
|
|---|
| 342 | d['energy range']=(h[57],h[58])
|
|---|
| 343 |
|
|---|
| 344 | d['kin. energy cutoff for hadrons in GeV'] = h[59]
|
|---|
| 345 | d['kin. energy cutoff for muons in GeV'] = h[60]
|
|---|
| 346 | d['kin. energy cutoff for electrons in GeV'] = h[61]
|
|---|
| 347 | d['energy cutoff for photons in GeV'] = h[62]
|
|---|
| 348 |
|
|---|
| 349 | d['NFLAIN'] = h[63]
|
|---|
| 350 | d['NFLDIF'] = h[64]
|
|---|
| 351 | d['NFLPI0'] = h[65]
|
|---|
| 352 | d['NFLPIF'] = h[66]
|
|---|
| 353 | d['NFLCHE'] = h[67]
|
|---|
| 354 | d['NFRAGM'] = h[68]
|
|---|
| 355 |
|
|---|
| 356 | d["Earth's magnetic field in uT: (x,z)"] = ( h[69], h[70] )
|
|---|
| 357 |
|
|---|
| 358 | d['flag for activating EGS4'] = h[71]
|
|---|
| 359 | d['flag for activating NKG'] = h[72]
|
|---|
| 360 | d['low-energy hadr. model flag (1.=GHEISHA, 2.=UrQMD, 3.=FLUKA)'] = h[73]
|
|---|
| 361 | d['high-energy hadr. model flag (0.=HDPM,1.=VENUS, 2.=SIBYLL,3.=QGSJET, 4.=DPMJET, 5.=NE X US, 6.=EPOS)'] = h[74]
|
|---|
| 362 | d['CERENKOV Flag (is a bitmap --> usersguide)'] = hex(int(round(h[75])))
|
|---|
| 363 | d['NEUTRINO flag'] = h[76]
|
|---|
| 364 | d['CURVED flag (0=standard, 2=CURVED)'] = h[77]
|
|---|
| 365 | d['computer flag (3=UNIX, 4=Macintosh)'] = h[78]
|
|---|
| 366 | d['theta interval (in degree): (lower, upper edge) '] = ( h[79], h[80] )
|
|---|
| 367 | d['phi interval (in degree): (lower, upper edge) '] = ( h[81], h[82] )
|
|---|
| 368 |
|
|---|
| 369 | d['Cherenkov bunch size in the case of Cherenkov calculations'] = h[83]
|
|---|
| 370 | d['number of Cherenkov detectors in (x, y) direction'] = (h[84], h[85])
|
|---|
| 371 | d['grid spacing of Cherenkov detectors in cm (x, y) direction'] = ( h[86], h[87])
|
|---|
| 372 | d['length of each Cherenkov detector in cm in (x, y) direction'] = ( h[88], h[89])
|
|---|
| 373 | d['Cherenkov output directed to particle output file (= 0.) or Cherenkov output file (= 1.)'] = h[90]
|
|---|
| 374 |
|
|---|
| 375 | d['angle (in rad) between array x-direction and magnetic north'] = h[91]
|
|---|
| 376 | d['flag for additional muon information on particle output file'] = h[92]
|
|---|
| 377 | d['step length factor for multiple scattering step length in EGS4'] = h[93]
|
|---|
| 378 | d['Cherenkov bandwidth in nm: (lower, upper) end'] = ( h[94], h[95] )
|
|---|
| 379 | d['number i of uses of each Cherenkov event'] = h[96]
|
|---|
| 380 |
|
|---|
| 381 | core_location_for_scattered_events_x_location = []
|
|---|
| 382 | core_location_for_scattered_events_y_location = []
|
|---|
| 383 | for i in range(97,117):
|
|---|
| 384 | core_location_for_scattered_events_x_location.append(h[i])
|
|---|
| 385 | for i in range(117,137):
|
|---|
| 386 | core_location_for_scattered_events_y_location.append(h[i])
|
|---|
| 387 | d['core location for scattered events in cm: (x,y)'] = zip(
|
|---|
| 388 | core_location_for_scattered_events_x_location,
|
|---|
| 389 | core_location_for_scattered_events_y_location)
|
|---|
| 390 |
|
|---|
| 391 | d['SIBYLL interaction flag (0.= no SIBYLL, 1.=vers.1.6; 2.=vers.2.1)'] = h[137]
|
|---|
| 392 | d['SIBYLL cross-section flag (0.= no SIBYLL, 1.=vers.1.6; 2.=vers.2.1)'] = h[138]
|
|---|
| 393 | d['QGSJET interact. flag (0.=no QGSJET, 1.=QGSJETOLD,2.=QGSJET01c, 3.=QGSJET-II)'] = h[139]
|
|---|
| 394 | d['QGSJET X-sect. flag (0.=no QGSJET, 1.=QGSJETOLD,2.=QGSJET01c, 3.=QGSJET-II)'] = h[140]
|
|---|
| 395 | d['DPMJET interaction flag (0.=no DPMJET, 1.=DPMJET)'] = h[141]
|
|---|
| 396 | d['DPMJET cross-section flag (0.=no DPMJET, 1.=DPMJET)'] = h[142]
|
|---|
| 397 | d['VENUS/NE X US/EPOS cross-section flag (0=neither, 1.=VENUSSIG,2./3.=NEXUSSIG, 4.=EPOSSIG)'] = h[143]
|
|---|
| 398 | d['muon multiple scattering flag (1.=Moliere, 0.=Gauss)'] = h[144]
|
|---|
| 399 | d['NKG radial distribution range in cm'] = h[145]
|
|---|
| 400 | d['EFRCTHN energy fraction of thinning level hadronic'] = h[146]
|
|---|
| 401 | d['EFRCTHN x THINRAT energy fraction of thinning level em-particles'] = h[147]
|
|---|
| 402 | d['actual weight limit WMAX for thinning hadronic'] = h[148]
|
|---|
| 403 | d['actual weight limit WMAX x WEITRAT for thinning em-particles'] = h[149]
|
|---|
| 404 | d['max. radius (in cm) for radial thinning'] = h[150]
|
|---|
| 405 | d['viewing cone VIEWCONE (in deg): (inner, outer) angle'] = (h[151], h[152])
|
|---|
| 406 | d['transition energy high-energy/low-energy model (in GeV)'] = h[153]
|
|---|
| 407 | d['skimming incidence flag (0.=standard, 1.=skimming)'] = h[154]
|
|---|
| 408 | d['altitude (cm) of horizontal shower axis (skimming incidence)'] = h[155]
|
|---|
| 409 | d['starting height (cm)'] = h[156]
|
|---|
| 410 | d['flag indicating that explicite charm generation is switched on'] = h[156]
|
|---|
| 411 | d['flag for hadron origin of electromagnetic subshower on particle tape'] = h[157]
|
|---|
| 412 | d['flag for observation level curvature (CURVOUT) (0.=flat, 1.=curved)'] = h[166]
|
|---|
| 413 |
|
|---|
| 414 | def printHeaderDict(self):
|
|---|
| 415 | print " Event Header "
|
|---|
| 416 | print "="*70
|
|---|
| 417 | pprint(self.header_dict)
|
|---|
| 418 | print "="*70
|
|---|
| 419 | print " End of: Event Header "
|
|---|
| 420 | print "="*70
|
|---|
| 421 |
|
|---|
| 422 | def __repr__(self):
|
|---|
| 423 | self.printInfoDict()
|
|---|
| 424 | return ""
|
|---|
| 425 |
|
|---|
| 426 |
|
|---|
| 427 | ##############################################################################
|
|---|
| 428 | ##############################################################################
|
|---|
| 429 | #
|
|---|
| 430 | # Two helper functions for plotting the contents of a corsika run follow
|
|---|
| 431 | #
|
|---|
| 432 |
|
|---|
| 433 | def plot_x_y(run, x, y, bins):
|
|---|
| 434 |
|
|---|
| 435 | plt.ion()
|
|---|
| 436 | fig = plt.figure()
|
|---|
| 437 |
|
|---|
| 438 | if y == None:
|
|---|
| 439 | # find min and max of x
|
|---|
| 440 | xmin = np.inf
|
|---|
| 441 | xmax = -np.inf
|
|---|
| 442 |
|
|---|
| 443 | for ev in run.events:
|
|---|
| 444 | if ev.photons[:,x].min() < xmin:
|
|---|
| 445 | xmin = ev.photons[:,x].min()
|
|---|
| 446 | if ev.photons[:,x].max() > xmax:
|
|---|
| 447 | xmax = ev.photons[:,x].max()
|
|---|
| 448 |
|
|---|
| 449 | ra = np.array([xmin,xmax])
|
|---|
| 450 | Hsum = None
|
|---|
| 451 | for ev in run.events:
|
|---|
| 452 | p = ev.photons
|
|---|
| 453 | H,xb = np.histogram(p[:,x], bins=bins, range=ra)
|
|---|
| 454 | if Hsum == None:
|
|---|
| 455 | Hsum = H
|
|---|
| 456 | else:
|
|---|
| 457 | Hsum += H
|
|---|
| 458 |
|
|---|
| 459 | width = 0.7*(xb[1]-xb[0])
|
|---|
| 460 | center = (xb[:-1]+xb[1:])/2
|
|---|
| 461 | plt.bar(center, Hsum, align = 'center', width = width)
|
|---|
| 462 |
|
|---|
| 463 | else:
|
|---|
| 464 |
|
|---|
| 465 |
|
|---|
| 466 | # find min and max of x and y
|
|---|
| 467 | xmin = np.inf
|
|---|
| 468 | xmax = -np.inf
|
|---|
| 469 | ymin = np.inf
|
|---|
| 470 | ymax = -np.inf
|
|---|
| 471 |
|
|---|
| 472 | for ev in run.events:
|
|---|
| 473 | #print x,".min=", ev.photons[:,x].min()
|
|---|
| 474 | #print x,".max=", ev.photons[:,x].max()
|
|---|
| 475 | #print y,".min=", ev.photons[:,y].min()
|
|---|
| 476 | #print y,".max=", ev.photons[:,y].max()
|
|---|
| 477 |
|
|---|
| 478 | if ev.photons is None: continue
|
|---|
| 479 |
|
|---|
| 480 |
|
|---|
| 481 | if ev.photons[:,x].min() < xmin:
|
|---|
| 482 | xmin = ev.photons[:,x].min()
|
|---|
| 483 | if ev.photons[:,x].max() > xmax:
|
|---|
| 484 | xmax = ev.photons[:,x].max()
|
|---|
| 485 | if ev.photons[:,y].min() < ymin:
|
|---|
| 486 | ymin = ev.photons[:,y].min()
|
|---|
| 487 | if ev.photons[:,y].max() > ymax:
|
|---|
| 488 | ymax = ev.photons[:,y].max()
|
|---|
| 489 |
|
|---|
| 490 | #xmin, ymin = -1, -1
|
|---|
| 491 | #xmax, ymax = 1,1
|
|---|
| 492 | # gather complete histo
|
|---|
| 493 | ra = np.array([xmin,xmax,ymin,ymax]).reshape(2,2)
|
|---|
| 494 | Hsum = None
|
|---|
| 495 | for ev in run.events:
|
|---|
| 496 | p = ev.photons
|
|---|
| 497 | if p is None: continue
|
|---|
| 498 | H,xb,yb = np.histogram2d(p[:,x], p[:,y],bins=[bins,bins], range=ra)
|
|---|
| 499 | if Hsum == None:
|
|---|
| 500 | Hsum = H
|
|---|
| 501 | else:
|
|---|
| 502 | Hsum += H
|
|---|
| 503 |
|
|---|
| 504 | extent = [ yb[0] , yb[-1] , xb[0] , xb[-1] ]
|
|---|
| 505 | plt.imshow(Hsum, extent=extent, interpolation='nearest')
|
|---|
| 506 |
|
|---|
| 507 | plt.colorbar()
|
|---|
| 508 |
|
|---|
| 509 |
|
|---|
| 510 | def plot_xy_rel_core_location(run, bins):
|
|---|
| 511 |
|
|---|
| 512 | plt.ion()
|
|---|
| 513 | fig = plt.figure()
|
|---|
| 514 |
|
|---|
| 515 |
|
|---|
| 516 | # find min and max of x and y
|
|---|
| 517 | xmin = np.inf
|
|---|
| 518 | xmax = -np.inf
|
|---|
| 519 | ymin = np.inf
|
|---|
| 520 | ymax = -np.inf
|
|---|
| 521 |
|
|---|
| 522 | for ev in run.events:
|
|---|
| 523 | #print x,".min=", ev.photons[:,x].min()
|
|---|
| 524 | #print x,".max=", ev.photons[:,x].max()
|
|---|
| 525 | #print y,".min=", ev.photons[:,y].min()
|
|---|
| 526 | #print y,".max=", ev.photons[:,y].max()
|
|---|
| 527 |
|
|---|
| 528 | if ev.photons is None: continue
|
|---|
| 529 |
|
|---|
| 530 | # get core-location from event header dictionary
|
|---|
| 531 | cl = np.array(ev.header_dict['core location for scattered events in cm: (x,y)'][0])
|
|---|
| 532 |
|
|---|
| 533 | xy = ev.photons[:,1:3] - cl
|
|---|
| 534 |
|
|---|
| 535 |
|
|---|
| 536 | if xy[:,0].min() < xmin:
|
|---|
| 537 | xmin = xy[:,0].min()
|
|---|
| 538 | if xy[:,0].max() > xmax:
|
|---|
| 539 | xmax = xy[:,0].max()
|
|---|
| 540 | if xy[:,1].min() < ymin:
|
|---|
| 541 | ymin = xy[:,1].min()
|
|---|
| 542 | if xy[:,1].max() > ymax:
|
|---|
| 543 | ymax = xy[:,1].max()
|
|---|
| 544 |
|
|---|
| 545 | #xmin, ymin = -1, -1
|
|---|
| 546 | #xmax, ymax = 1,1
|
|---|
| 547 | # gather complete histo
|
|---|
| 548 | ra = np.array([xmin,xmax,ymin,ymax]).reshape(2,2)
|
|---|
| 549 | Hsum = None
|
|---|
| 550 | for index, ev in enumerate(run.events):
|
|---|
| 551 | if ev.photons is None: continue
|
|---|
| 552 | # get core-location from event header dictionary
|
|---|
| 553 | cl = np.array(ev.header_dict['core location for scattered events in cm: (x,y)'][0])
|
|---|
| 554 | xy = ev.photons[:,1:3] - cl
|
|---|
| 555 |
|
|---|
| 556 | H,xb,yb = np.histogram2d(xy[:,0], xy[:,1], bins=[bins,bins], range=ra)
|
|---|
| 557 | if Hsum == None:
|
|---|
| 558 | Hsum = H
|
|---|
| 559 | else:
|
|---|
| 560 | Hsum += H
|
|---|
| 561 | #print "index: ", index, "H.sum:", H.sum()
|
|---|
| 562 | #print "index: ", index, "Hsum.sum:", Hsum.sum()
|
|---|
| 563 |
|
|---|
| 564 | extent = [ yb[0] , yb[-1] , xb[0] , xb[-1] ]
|
|---|
| 565 | plt.imshow(Hsum, extent=extent, interpolation='nearest')
|
|---|
| 566 | plt.colorbar()
|
|---|
| 567 | return Hsum
|
|---|
| 568 |
|
|---|
| 569 |
|
|---|
| 570 |
|
|---|
| 571 |
|
|---|
| 572 |
|
|---|
| 573 |
|
|---|
| 574 |
|
|---|
| 575 |
|
|---|
| 576 | def read_corsika_file( filename ):
|
|---|
| 577 | f = open(filename, 'rb')
|
|---|
| 578 |
|
|---|
| 579 | while True: # the only break out of this while is finding a run footer block - "RUNE"
|
|---|
| 580 |
|
|---|
| 581 | # Each "block" in a corika file is always 273 "words" long.
|
|---|
| 582 | # a work is a 4-byte value, usually a double, but might also be an int.
|
|---|
| 583 | # The first "word" of a block is sometimes a sequence of 4 ascii
|
|---|
| 584 | # charaters, in order to tell the reader, what type of block follows next.
|
|---|
| 585 | # these sequences can be:
|
|---|
| 586 | # * "RUNH" - a run header block
|
|---|
| 587 | # * "RUNE" - a run footer block
|
|---|
| 588 | # * "EVTH" - an event header block
|
|---|
| 589 | # * "EVTE" - an event footer block
|
|---|
| 590 | # the corsika manual states more of there, but those are not
|
|---|
| 591 | # expected in a corsika "cer-file".
|
|---|
| 592 | block = f.read(273*4)
|
|---|
| 593 | if not block:
|
|---|
| 594 | raise Exception('End of file '+ filename +' reached, before ' +
|
|---|
| 595 | '"RUNE" word found. File might be corrupted.')
|
|---|
| 596 | if not len(block) == 273*4:
|
|---|
| 597 | raise Exception("""While reading %s a block smaller than 273*4 bytes
|
|---|
| 598 | was read, this might be due to non-blocking read
|
|---|
| 599 | which is clearly a bug in this software.
|
|---|
| 600 | you might want to simply try again. """ % (filename,))
|
|---|
| 601 |
|
|---|
| 602 | # what ever block we have, we treat the first 4 bytes as
|
|---|
| 603 | # a 4-character string. In case it is not one of the
|
|---|
| 604 | # "block descriptors", we assume it is a data-block.
|
|---|
| 605 | block_type = block[0:4]
|
|---|
| 606 |
|
|---|
| 607 | # this reader has no error handling of the underlying file at all.
|
|---|
| 608 | # so in case there is an EVTH found before a RUNH the reader will
|
|---|
| 609 | # probably totally freak out :-)
|
|---|
| 610 | # This is a task for captain state-machine :-)
|
|---|
| 611 |
|
|---|
| 612 | if block_type == 'RUNH':
|
|---|
| 613 | # if we hava a run-header block, the trailing 272 words are
|
|---|
| 614 | # unpacked into a list of floats. from this list a Run-object
|
|---|
| 615 | # is created, which is mainly empty up to that point.
|
|---|
| 616 | run_header = struct.unpack( '272f', block[4:] )
|
|---|
| 617 | run = Run(run_header)
|
|---|
| 618 |
|
|---|
| 619 | elif block_type == 'EVTH':
|
|---|
| 620 | # directly after the RUNH an event header is expected.
|
|---|
| 621 | # so if that event header is found. again the trailing
|
|---|
| 622 | # 272 words areunpacked into a list of 272 floats and used to
|
|---|
| 623 | # call the constructor of an Event instance.
|
|---|
| 624 | # this event, still it is mainly empty, is appended to the
|
|---|
| 625 | # list of events, maintained by the run object.
|
|---|
| 626 | #
|
|---|
| 627 | # go on reading in the else-block of this of condition.
|
|---|
| 628 | event_header = struct.unpack( '272f', block[4:] )
|
|---|
| 629 | event = Event(event_header)
|
|---|
| 630 | run.events.append(event)
|
|---|
| 631 |
|
|---|
| 632 | elif block_type == 'EVTE':
|
|---|
| 633 | # after all data blocks for a certain event have been found
|
|---|
| 634 | # we will eventually find a event footer block.
|
|---|
| 635 | # the footer does not contain much information, but it is still
|
|---|
| 636 | # safed with the current event.
|
|---|
| 637 | # It could be used to check, if we found a footer, which fits to the
|
|---|
| 638 | # event header, but that is not done here.
|
|---|
| 639 | event.footer = struct.unpack( '272f', block[4:] )
|
|---|
| 640 |
|
|---|
| 641 | # When footer is found, we know that the event is
|
|---|
| 642 | # fully read. So we can make a python dict, which contains
|
|---|
| 643 | # some short info about his event, like the number of photons
|
|---|
| 644 | # which is not stored in the event header.
|
|---|
| 645 | event.makeInfoDict()
|
|---|
| 646 | event = None
|
|---|
| 647 |
|
|---|
| 648 | elif block_type == 'RUNE':
|
|---|
| 649 | # when the run footer block is found, we know we have
|
|---|
| 650 | # reached the end of the file, since as far as I have seen
|
|---|
| 651 | # there is only one run in one corsika file.
|
|---|
| 652 | # In addition one can find a lot of zeroes at the end of a corsika
|
|---|
| 653 | # file after the run footer block, so it is not safe to
|
|---|
| 654 | # simply read until the end of file and try to parse everything.
|
|---|
| 655 | run.footer = struct.unpack( '272f', block[4:] )
|
|---|
| 656 | break
|
|---|
| 657 | else:
|
|---|
| 658 | # since there was no block-descripter found, we assume this
|
|---|
| 659 | # block is a data block. the 273 words of this block are
|
|---|
| 660 | # again unpacked into 273 floats this time.
|
|---|
| 661 | # from the corsika manual we know, that a data block is grouped
|
|---|
| 662 | # into 39 groups of 7 words, so we form an np.array of shape
|
|---|
| 663 | # 39 x 7
|
|---|
| 664 | data_block = struct.unpack( '273f', block )
|
|---|
| 665 | data_array = np.array( data_block, dtype=np.float32 ).reshape(39,7)
|
|---|
| 666 |
|
|---|
| 667 | # the data we read in here, does not comply with the original
|
|---|
| 668 | # data sheet.
|
|---|
| 669 | # corsika was edited by somebody calles Dorota from Lodz in order to write the
|
|---|
| 670 | # wavelength of the simulated cherenkov photons and some other info.
|
|---|
| 671 | # I did not find any docu, but had a look into the code, which
|
|---|
| 672 | # luckily is stored right next to the binary. :-)
|
|---|
| 673 | # according to that code the 1st element of the 7-element-tuple
|
|---|
| 674 | # belonging to each cherenkov-photon stores 3 types of information:
|
|---|
| 675 | # a parameter called J ... which might be the parent particles type
|
|---|
| 676 | # a parameter called IMOV, which seems to be constant = 1
|
|---|
| 677 | # and which seems to have something to do with reusing showers
|
|---|
| 678 | # a parameter called WAVELENGTH which seems to be the cherekov
|
|---|
| 679 | # photon wavelength in nm.
|
|---|
| 680 | #
|
|---|
| 681 | # the code looks like this:
|
|---|
| 682 | # J*100000. + IMOV*1000. + WAVELENGTH
|
|---|
| 683 | new_extra_data = np.zeros( (39,3), dtype=np.float32 )
|
|---|
| 684 | for index, row in enumerate(data_array):
|
|---|
| 685 | integer = long(round(row[0]))
|
|---|
| 686 | J = integer/100000
|
|---|
| 687 | IMOV = (integer-J*100000)/1000
|
|---|
| 688 | WAVELENGTH = np.mod(row[0], 1000.)
|
|---|
| 689 | new_extra_data[index,0] = J
|
|---|
| 690 | new_extra_data[index,1] = IMOV
|
|---|
| 691 | new_extra_data[index,2] = WAVELENGTH
|
|---|
| 692 |
|
|---|
| 693 | data_array = np.concatenate( (data_array,new_extra_data), axis=1)
|
|---|
| 694 |
|
|---|
| 695 | # so at this point we have made from the 39x7 shaped array an
|
|---|
| 696 | # 39x10 shaped array.
|
|---|
| 697 | # the contents of this array are 4 byte floats and have
|
|---|
| 698 | # the following meansings:
|
|---|
| 699 | # d = data_array[0]
|
|---|
| 700 | # d[0] - encoded info from dorota. see index:7,8,9
|
|---|
| 701 | # d[1] - x position in cm; x-axis points north in corsika
|
|---|
| 702 | # d[2] - y position in cm; y-axis points west in corsika
|
|---|
| 703 | # both positions are measured at the "observation level" aka ground
|
|---|
| 704 | # d[3] - direction cosine to x-axis, called "u" in corsika
|
|---|
| 705 | # d[4] - direction cosine to y-axis, called "v" in corsika
|
|---|
| 706 | # d[5] - time since first interaction in ns
|
|---|
| 707 | # d[6] - height of production in cm;
|
|---|
| 708 | #
|
|---|
| 709 | # now we can't be sure, since we have it reverse engineered.
|
|---|
| 710 | # d[7] - "J", maybe particle ID of ch. light emitting particle
|
|---|
| 711 | # d[8] - "IMOV", apparently always 1, connection with reusing of showers
|
|---|
| 712 | # d[9] - "WAVELENGTH" wavelength of the photon in nm.
|
|---|
| 713 |
|
|---|
| 714 |
|
|---|
| 715 | # since a data block has a fixed size, it can contain
|
|---|
| 716 | # empty rows, we want to find the first empty row here.
|
|---|
| 717 | first_empty_row_index = False
|
|---|
| 718 | for row_index, row in enumerate(data_array):
|
|---|
| 719 | if (row == 0.).all():
|
|---|
| 720 | first_empty_row_index = row_index
|
|---|
| 721 | break
|
|---|
| 722 |
|
|---|
| 723 | # and here we simply cut off all empty rows.
|
|---|
| 724 | if first_empty_row_index:
|
|---|
| 725 | data_array = data_array[ 0:first_empty_row_index, : ]
|
|---|
| 726 |
|
|---|
| 727 | # here we append our data_array to the np.array, which is stored
|
|---|
| 728 | # in the current event. keep in mind that the current event can simply
|
|---|
| 729 | # be retrieved by "event".
|
|---|
| 730 | if event.photons == None:
|
|---|
| 731 | event.photons = data_array
|
|---|
| 732 | else:
|
|---|
| 733 | event.photons = np.concatenate( (event.photons, data_array), axis=0 )
|
|---|
| 734 |
|
|---|
| 735 | f.close()
|
|---|
| 736 |
|
|---|
| 737 | run = add_all_photon_table(run)
|
|---|
| 738 |
|
|---|
| 739 | return run
|
|---|
| 740 |
|
|---|
| 741 |
|
|---|
| 742 |
|
|---|
| 743 |
|
|---|
| 744 | if __name__ == '__main__':
|
|---|
| 745 |
|
|---|
| 746 | run = read_corsika_file( sys.argv[1] )
|
|---|
| 747 | print 'Run object "run" created.'
|
|---|
| 748 | print 'type: print run'
|
|---|
| 749 | print 'or'
|
|---|
| 750 | print 'for event in run.events:'
|
|---|
| 751 | print ' for photon in event.photons:'
|
|---|
| 752 | print ' print photon'
|
|---|
| 753 | print
|
|---|
| 754 | print 'to loop over everyting...'
|
|---|
| 755 |
|
|---|
| 756 |
|
|---|
| 757 |
|
|---|
| 758 | #plot_x_y(run, 3,4,1000)
|
|---|
| 759 | #plt.title("photon distribution ... in theta vs phi ... or other way round")
|
|---|
| 760 | #plot_x_y(run, 1,2,1000)
|
|---|
| 761 | #plt.title("photon distribution ... x vs y .. or other way round")
|
|---|
| 762 |
|
|---|
| 763 | #Hsum = plot_xy_rel_core_location(run, 100)
|
|---|
| 764 | #plt.title("photon distribution 100 bins; file:" + sys.argv[1])
|
|---|
| 765 | #plt.xlabel('photon y location relative to "core location" in cm')
|
|---|
| 766 | #plt.ylabel('photon x location relative to "core location" in cm')
|
|---|
| 767 |
|
|---|
