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