Corsika Parse Output
About
本文将会是一个简单的踩坑介绍.
Corsika Install and Run
发邮件, 下载程序, 使用 ./coconut 进行安装, 最后 run 目录下开始运行.
Corsika Output
Corsika 程序的输出 DAT<run-number> 是一个 Fortran 的二进制结构.
本来还有一个比较友好的数据读取和处理程序 COAST.
但是一个比较尴尬的事情是貌似它有些问题无法正常读取数据 (学长说的),
在我的电脑上无法正常编译 (configure 文件没了, make 爆了).
那么与其去修代码解决依赖问题, 我决定听从学长的建议, 写一个程序自己读.
找到官方的文档 (CORSIKA_GUIDE), 然后找到 Output 部分. 嗯, 接下来就是一些肮脏的过程了.
也不算坑的一些东西
Corsika 的安装和使用
二进制数据读取
我用的是 lisp-binary 来帮助我从二进制文件里面读取数据. 毕竟不是严格意义上搞这些的人, 啥 IEEE 的浮点数定义和实现也不是很懂. 总之能读就好了吧.
按照 Corsika Guide 里面定义的块结构, 可以如下定义如何读取数据:
(defbinary particle-data (:export t :byte-order :little-endian)
(description 0d0 :type single-float)
(p-x 0d0 :type single-float)
(p-y 0d0 :type single-float)
(p-z 0d0 :type single-float)
(x 0d0 :type single-float)
(y 0d0 :type single-float)
(time 0d0 :type single-float))
(defbinary particle-data-block (:export t :byte-order :little-endian)
(particles #() :type (simple-array particle-data (38))))
然后通过 (read-binary 'particle-data-block stream) 从输入流中读取信息.
当然, 这里有一些小小的 Trick, 在之后会介绍.
Sub Block Read and Unread
一个简单的读取想法就是读数, 判断类型, 然后根据类型返回结构. 道理是这个样子, 但是有一个小小问题, 因为 lisp-binary 是用 stream 的方式进行数据读取的, 读完了之后, 已经读过的东西是不能再读的.
除非你把它放到一个缓冲块里面. (注: 我没有学过更加高深的数据读取方法, 所以我也说不好这样的操作是否最优, 至少当前我觉得能用就好了)
于是我写了这样的一个宏来帮助我实现 unread 这样的操作:
(defmacro with-stream-block ((stream block) &body body)
(alexandria:with-gensyms (in out)
`(flexi-streams:with-input-from-sequence
(,in (flexi-streams:with-output-to-sequence (,out)
(write-binary ,block ,out)))
(with-wrapped-in-bit-stream (,stream ,in)
,@body))))
(defun read-block (corsika)
(let* ((dummy (read-binary 'dummy-sub-block corsika))
(type (dummy-sub-block-type dummy)))
(values
(with-stream-block (stream dummy)
(cond ((string= "RUNH" type) (read-binary 'run-header stream))
((string= "EVTH" type) (read-binary 'event-header stream))
((string= "LONG" type) (read-binary 'longitudinal-block stream))
((string= "EVTE" type) (read-binary 'event-end stream))
((string= "RUNE" type) (read-binary 'run-end stream))
(T (read-binary 'particle-data-block stream))))
type)))
这里默认 DATA BLOCK 全部都是粒子数据块.
并且会发现这里是一口气读取整个 sub-block, 并且只使用头部一个很小的部分进行类型判断, 这样有点慢. 最好是能够根据头部信息进行动态判断, 但是这个 particle-data-block 的结构, 和其他的 sub-block 的结构并不一样, 所以还是有点坑.
Data Block
首先到网上找了一个讲解 Corsika 的教学 PPT (Corsika_Installation_Physics.pdf), 然后里面的一节介绍了 Corsika 的二进制数据结构:
(当然, 更加详细的可以看 Corsika Guide 在第 10 节的 Output 里面的说明)
对于 THIN 参数, 可以跳过不理会先 (我应该没有在安装的过程中选择这个模式).
每个块的大小和作用都讲得挺清楚的, 所以看看图应该就能够理解.
这里有一个让人一开始比较容易误会的地方就是:
1 block = 5733(6552) words (4 bytes) = 21 sub-blocks of 273(312) words
并且在 Corsika Guide 里面也是这样画的:
看起来一个 Block 的结构里面就是类似于:
Run Header
Event Header
Data Block
...
Long Block
Event End
...
Run End
这样的东西, 那么如果真就这样去读取数据的话, 那么就进入了一个比较坑爹的大坑了.
尤其是当你用二进制读取程序去看这个最后输出的文件的话,
发现好像这个 RUN HEADER 和 RUN END 只有一对 (对于一次的 RUN),
但是不同事件数 EVENT 却有不同的长度, 那么这个 1 block = 21 sub-blocks 究竟是啥?
中间我试过去掉开头的 32 bytes (学长建议), 后来发现开头的 32 bytes 在无符号整数下代表
block 的大小 (对于我手上的, 就是 22932), 紧跟着的就是子块结构, 所以我的做法就是:
;;; This is not TRUE!!!!
(defun parse-corsika (corsika)
;; how to read block shall depend on output type
(read-corsika-type corsika)
(loop for sub-block = (read-block corsika)
while (not (eq 'run-end (type-of sub-block)))
collect sub-block))
这样读到的结果是什么呢? 答案是前几个可以正常地读出来, 而之后的几个就开始出现乱码了. (其实不是乱码, 而是输出很明显不符合物理), 并且一只读下去会卡死并报出文件空的读取错误. 显然事情并没有那么简单.
那么最终的答案是什么呢?
答案是, 这个 block 就是个没有任何意义, 也没有啥结构的, 强行对输出文件每 21 个 sub-block 进行一划分的一个东西. (可能来源于历史上的卡片输出, 但是对于现在的这个, 真的是, 没啥必要). 也就是真实的数据如下:
block marker (22932)
RUN HEADER
EVENT HEADER
DATA BLOCK
... (21 subblocks)
block marker (22932)
block marker (22932)
DATA BLOCK
LONG BLOCK
...
EVENT END
... (21 sub-blocks)
block marker (22932)
block marker (22932)
...
RUN END
... (empty with 0)
block marker (22932)
如果你用下面的程序来读取的话:
(with-open-binary-file (corsika output-data-path)
(apply #'append
(loop for i from 0 below 5
do (read-corsika-type corsika)
collect (loop for j below 21
for sub-block = (read-block corsika)
collect sub-block)
do (read-corsika-type corsika))))
那么就能合理的看到最终的输出了
(#S(RUN-HEADER
:TYPE "RUNH"
:RUN-NUMBER 69.0
:START-DATE 231009.0
:PROGRAM-VERSION 7.75
;;; 略
:NFLCHE+100NFRAGM 200.0)
#S(EVENT-HEADER
:TYPE "EVTH"
:EVENT-NUMBER 1.0
:PARTICLE-ID 5626.0
;;; 略
:NO-USE #(0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0))
#S(PARTICLE-DATA-BLOCK
:PARTICLES #(#S(PARTICLE-DATA
:DESCRIPTION 6561.0
:P-X -2.552703
:P-Y 0.38220033
:P-Z 1.9585683
:X -152077.31
:Y 41590.297
:TIME 176972.77)
;;; 略
)))
Parse Data Structure
一个简单的数据读取程序如下:
(defun parse-corsika (file)
(with-open-binary-file (corsika file)
(let ((eof NIL))
(apply #'nconc
(loop do (read-corsika-type corsika)
collect (loop for counter below 21
for sub-block = (read-block corsika)
collect sub-block
while (not (eq 'run-end (type-of sub-block)))
finally (if (eq 'run-end (type-of sub-block))
(setf eof T)))
while (not eof)
do (read-corsika-type corsika))))))
缺点就是有点太慢了, 并且出来的东西还不是有层次结构的东西, 需要能够更快一点, 更好一点…
写了一个简单的 parser, 看看能不能读出结构:
(defun parse-corsika-structure (corsika)
(parse-run-block corsika))
(defun parse-run-block (corsika)
(let ((header (car corsika))
(remain (cdr corsika)))
(when (typep header 'run-header)
(list :header header
:events (loop for event = (multiple-value-bind (event rest)
(parse-event-block remain)
(setf remain rest)
event)
while event
collect event)
:end (car remain)))))
(defun parse-event-block (corsika)
(let ((header (car corsika))
(remain (cdr corsika)))
(if (typep header 'event-header)
(values
(list :header header
:data-block (loop for data-block = (car remain)
while (typep data-block 'particle-data-block)
do (setf remain (cdr remain))
collect data-block)
:long-block (loop for long-block = (car remain)
while (typep long-block 'longitudinal-block)
do (setf remain (cdr remain))
collect long-block)
:end (car remain))
(cdr remain))
(values NIL corsika))))
用一个简单的函数来处理一下最终的结果, 防止出来的数据太多把屏幕给挤爆:
(defun map-tree (function tree)
(if (listp tree)
(loop for elem in tree
collect (map-tree function elem))
(funcall function tree)))
(map-tree (lambda (elem)
(if (keywordp elem) elem (type-of elem)))
(parse-corsika-structure
(parse-corsika file)))
结果速览
(:HEADER RUN-HEADER :EVENTS
((:HEADER EVENT-HEADER :DATA-BLOCK
(PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK
PARTICLE-DATA-BLOCK PARTICLE-DATA-BLOCK)
:LONG-BLOCK
(LONGITUDINAL-BLOCK LONGITUDINAL-BLOCK LONGITUDINAL-BLOCK
LONGITUDINAL-BLOCK)
:END EVENT-END))
:END RUN-END)
又: 好像也没有那么慢, 因为之前读的是一个产生了 1000 个 EVENT 的模拟结果, 所以会读得慢一点, 对于一个 EVENT 的结果用时大约是 0.117 秒, 还算可以接受, 那么就不管了. 对于 19 个 EVENT 的结果, 大约是 0.647 秒, 也不是不行.
整理一下进入下一个阶段.
End
总之目前先这样, 之后再想想看有没有更加好的方法来解决这个问题吧… 还是先着重处理物理上的问题.