宇宙线的能量跨越10多个量级，流量强度跨越30多个量级，能谱近似符合dN/dE∝E-γ，其中，谱指数γ的取值约为2.7[1]。MeV~400 TeV能区宇宙线流量强度较大，可以用卫星直接探测，更高能区宇宙线由于流量强度太低，主要利用地面实验间接探测[1-2]。水切伦科夫探测器是地面宇宙线探测技术的一种，它具有全天候、宽视场等优势，主要通过测量原初宇宙线进入大气后产生的广延大气簇射重建原初宇宙线信息[1, 3]。位于墨西哥的高海拔水切伦科夫天文台(the High Altitude Water Cherenkov observatory, HAWC) 和位于四川稻城的高海拔宇宙线观测站-水切伦科夫探测器阵列(Large High Altitude Air Shower Observatory-Water Cherenkov Detector Array, LHAASO-WCDA) 等多个宇宙线探测实验都采用了水切伦科夫探测器技术。
广延大气簇射[4-5]是高能原初宇宙线进入大气层后与大气成份发生相互作用，通过强子级联和电磁级联产生大量次级粒子的过程，通常用纵向发展和横向分布来描述其物理特征。常用于描述纵向发展的模型有Gaisser-Hillas function[6-7]，Greisen function[7-8]和Gaussian-in-Age[7, 9]，三者都与实验测量值相吻合。广延大气簇射的横向分布通常用横向分布函数(Lateral Distribution Function, LDF)和NKG(Nishimura-Kamata-Greisen) 函数描述[4, 7, 10-11]。LDF是简单导出的、纯数学拟合的、描述平均分布的经验公式，它的参数不直接关联广延大气簇射的物理特性[5, 11-12]。NKG函数比LDF更接近物理本质，一方面NKG函数能够反映广延大气簇射的横向分布和纵向发展等物理特征，另一方面在拟合过程中还充分考虑了广延大气簇射的统计涨落[4-5, 7, 10]。
地面宇宙线探测实验通过探测包括电子(e+、e-)、缪子(μ-、μ+) 和光子在内的广延大气簇射次级粒子密度分布和前锋面到达时间反演事例的芯位、方向和能量等信息[10, 13-14]，这个过程称为宇宙线事例重建或广延大气簇射重建。芯位重建是广延大气簇射重建的重要内容，常用的芯位重建方法有重心法、树形分析法、最大似然法、二维高斯拟合法等，这些方法各自有独特的优势和不足[4, 15-17]。重心法和树形分析法的优势在于计算量小，节省计算资源；不足在于重建时破坏了事例的对称性特征。芯位重建的精度依赖重建出的芯位到探测器中心的距离。最大似然法重建芯位精度高，但是精度也依赖探测器单元触发数量(Nhit)，当Nhit＞1 000时芯位分辨率＜2 m，
精度相比前两者显著提高；但是一方面该方法计算量较大，需要占用大量计算资源，另一方面其精度依赖于能量、方向和探测器单元触发数量等诸多因素。二维高斯拟合法是一种简单的数学拟合方法，克服了重心法和树状分析法的部分缺点，对于芯位位于探测器阵列内的事例(芯内事例) 重建效果较好，而对于芯位位于探测器阵列外的事例不能准确重建[4, 17-19]。除最大似然法外，这些重建方法都需要通过筛选条件丢弃大量数据，只保留芯位在阵列中心附近小范围的事例，以牺牲探测器有效面积的方式来保证重建结果的可靠性，造成大量观测数据的浪费[10, 15]。此外水切伦科夫探测器数据的低信噪比也是重建的挑战之一，如LHAASO-WCDA的观测数据在去噪之前信噪比约为70/218，通过遍历法去噪之后信噪比约为22/70。也就是说对于触发70个探测器单元的事例，其中约有22个通道的信号是噪声，这些噪声给精确重建带来很大的困难，而上述方法都不能很好地消除噪声带来的影响[15]。
本文基于广延大气簇射的对称性和水切伦科夫探测器阵列的物理特征，提出了一种新的芯位重建方法——椭圆拟合法。该方法充分利用广延大气簇射的对称性，消除噪声对芯位重建的影响。在满足芯位重建精度的前提下，椭圆拟合法能够重建芯位落在探测器边沿的事例和芯外事例，提高探测器的事例利用率。
1 水切伦科夫探测器技术水切伦科夫探测器利用纯水作为介质，通过光电倍增管收集带电粒子在水中运动产生的切伦科夫光子来测量广延大气簇射。现有实验采用的水切伦科夫探测器阵列的单元探测器构造和探测原理基本相同，探测器单元的各种性能参数和布局存在一定的差别。LHAASO-WCDA是最具代表性的水切伦科夫探测器阵列之一。当宇宙线粒子在视场内产生广延大气簇射，次级粒子进入阵列内的纯水中做超光速(介质中的光速) 运动激发切伦科夫光，每个基本单元的光电倍增管收集的切伦科夫光电子数(NPE) 与该位置的广延大气簇射次级粒子密度成正比。探测器单元触发时间可以精确到纳秒量级，通过这些数据可以重建芯位、方向等信息[1, 20-21]。其他实验的水切伦科夫探测器阵列如HAWC的结构和布局与LHAASO-WCDA有所不同，但是探测原理基本一致[3]。LHAASO-WCDA由3个水池组成，具有3 120个单元探测器以及后级的电子学、定时、数据采集、触发判选、数据处理、标定等功能系统。单元探测器(图 1) 为5 m × 5 m的水域，深度为4.4 m，两个单元之间用隔光帘隔开, 避免来自同一次级粒子，尤其是缪子信号的串扰。WCDA采用3 120支大尺寸光电倍增管分别布设于每个单元的中央，置于水底，向上观测。另外我们还在水池中放置3 120支小尺寸光电倍增管，小尺寸光电倍增管放置在大尺寸光电倍增管旁边，用于扩大簇射粒子数测量的动态范围，从而实现高能宇宙线的高精度测量[1, 20-21]。本文以150 m × 150 m的LHAASO-WCDA1水切伦科夫探测器阵列为样本，阐述椭圆拟合法重建芯位的物理思路和算法设计。
|
| 图 1 (a) LHAASO-WCDA单元探测器结构示意图；(b) 探测示意图 Fig. 1 (a) Schematic diagram of LHAASO-WCDA unit detector structure; (b) detection diagram |
如图 1 (b)，当广延大气簇射以一定的天顶角入射到地面，到达探测器阵列平面，广延大气簇射次级粒子集中在很薄的前锋面内，次级粒子横向分布可以用NKG函数描述[4-5, 7, 10, 12-13, 22]：
| $ \rho (r) = \frac{{{N_{\rm{e}}}}}{{r_{\rm{M}}^2}}\frac{{\mathit{\Gamma }(4.5 - s)}}{{2\pi \mathit{\Gamma }(s)\mathit{\Gamma }(4.5 - 2s)}}{\left( {\frac{{{r^\prime }}}{{{r_{\rm{M}}}}}} \right)^{s - 2}}{\left( {1 + \frac{r}{{{r_{\rm{M}}}}}} \right)^{s - 4.5}}, $ | (1) |
其中，Ne为广延大气簇射的总粒子数，主要依赖于原初宇宙线粒子的能量；rM为莫里哀半径，与原初宇宙线的成分有关；s为描述广延大气簇射发展阶段的参数，称为广延大气簇射的年龄，当s=1时，广延大气簇射发展到极大
| $ \sigma(x, y)=\rho(r) \cos \varphi, $ | (2) |
其中，
| $ r=\frac{\left(x-x_0\right)^2\left(1-\sin ^2 \varphi \cos ^2 \theta\right)+\left(y-y_0\right)^2\left(1-\sin ^2 \varphi \sin ^2 \theta\right)-2\left(x-x_0\right)\left(y-y_0\right) \sin ^2 \varphi \sin \theta \cos \theta}{\sqrt{\left(x-x_0\right)^2+\left(y-y_0\right)^2}} ; $ | (3) |
(x0, y0) 为地平面上的芯位；(θ, φ) 为广延大气簇射的方向，θ为方位角，φ为天顶角；(x, y) 为平面上任意位置探测器单元的坐标。从表达式可以看出，σ(x, y) 与广延大气簇射的芯位和方向有关，其数学结构十分复杂，不方便分析其特征。为了更加直观地理解前锋面次级粒子的分布情况，我们进行了数值计算，并得到不同天顶角和方位角事例投射到地面的次级粒子分布情况，如图 2。
|
| 图 2 利用NKG函数数值计算得到的不同天顶角广延大气簇射带电次级粒子密度，方位角都是45°，其中, (a), (b)和(c) 的天顶角分别为0°，30°和45°，色度图为次级粒子密度的常用对数，(d), (e)和(f)分别为(a), (b)和(c) 的等高线表示 Fig. 2 The azimuth angles of the charged secondary particle density of EAS at different zenith angles calculated by using the NKG function are all 45°, where the zenith angles of (a), (b), and (c) are 0°, 30°, and 45°, respectively. The degree diagram is the common logarithm of the secondary particle density, and (d), (e), and (f) are represented by the contour lines of (a), (b), and (c), respectively |
从图 2可以看出，广延大气簇射次级粒子密度呈现明显的对称性，当天顶角为0°时次级粒子密度等高线为同心圆，当天顶角不为0时等高线为同心椭圆，天顶角越大，椭圆的离心率越大。因此，可以用椭圆系方程描述次级粒子密度等高线：
| $ \frac{\left[\left(x_n-x_0\right) \cos \theta+\left(y_n-y_0\right) \sin \theta\right]^2}{a_n{ }^2}+\frac{\left[\left(x_n-x_0\right) \sin \theta+\left(y_n-y_0\right) \cos \theta\right]^2}{b_n{ }^2}=1, $ | (4) |
其中，(x0, y0) 为一系列椭圆的共同中心；an和bn分别为第n个椭圆方程的长轴和短轴；θ为广延大气簇射的方位角；椭圆长短轴之比an/bn=f(φ, s)，是广延大气簇射天顶角φ和年龄s的函数。
水切伦科夫探测器阵列各基本单元光电倍增管探测的光电子数与该位置的次级粒子密度之间是线性关系，因此，光电子数等高线可以近似为椭圆。用椭圆方程拟合光电子数的等高线可以得到一系列的同心椭圆，这些椭圆的共同中心就是广延大气簇射的芯位，椭圆的长短轴之比反映广延大气簇射天顶角信息，椭圆的长轴反映广延大气簇射的方位角信息。拟合光电子数的等高线系列来重建芯位，一方面对于芯位位于探测器边沿使得广延大气簇射对称性严重破坏的事例，可以通过局部信息来“补齐”椭圆等高线(如图 3)，利用对称性重建准确的芯位；另一方面利用所有椭圆应该“同心”的特点，舍弃信噪比较低的数据，从而达到抑制噪声的效果。具体算法流程如图 4。
|
| 图 3 椭圆拟合法重建芯位示意图 Fig. 3 Schematic diagram of core position reconstruction by ellipse fitting method |
|
| 图 4 椭圆拟合法算法流程图 Fig. 4 Algorithm flow chart of ellipse fitting method |
(1) 读取事例获取对芯位重建有用的单元探测器坐标(xi, yi)；
(2) 把数据按照NPE降序排列，并把数据分成n组；
(3) 拟合得到n个椭圆方程，得到n个(x0i, y0i)；
(4) 求出其平均值x0, y0，标准差x0SD，y0SD和偏差值|x0-x0i|, |y0-y0i|；
(5) 判断是否满足|x0-x0i| < decut & |y0-y0i| < decut，如果满足进入下一步，否则去掉偏差值最大的点，回到第4步，直到满足条件或者n < 3；
(6) 计算芯位30 m内的圆形区域与30~40 m环形区域的NPE总和之比NBN，并且判断NBN是否大于NBNmin，若大于则输出上一步得到的(x0, y0)为重建芯位，否则判断重建失败，该事例不纳入统计。
该方法涉及的几个参数的选择和使用：
(1) 光电子数最大的点。对于芯内事例，光电子数最大的点不能用于拟合椭圆。这个点很可能位于椭圆中心附近。
(2) 椭圆系中拟合椭圆所用的数据点的个数。离芯位越远的椭圆等高线上的探测器单元越多，因此用于拟合一个椭圆的点越多，并且理论上与离芯位距离成正比。因此，我们采用从内到外线性增加的分组方式对数据进行分组。上文的推导表明拟合平面上任意椭圆最少需要5个点，故取每组数据为5n(n=1, 2, 3…)，即第1组数据5个点，第2组数据10个点，第3组15个点，以此类推。由此可得触发探测器单元数Nhit与拟合椭圆数n的关系为
| $ \frac{5 n(n+1)}{2} <N_{\text {hit }} \text {. } $ | (5) |
(3) decut。椭圆拟合法的核心是用椭圆方程拟合光电子数等高线得到椭圆中心(x0, y0) 作为宇宙线事例的芯位。由于所有的椭圆等高线应该是同心的，理论上只需要拟合一条等高线就可以得到芯位，但是为了减小统计误差，提高芯位重建的精度，我们需要充分利用数据拟合尽可能多的椭圆。由于物理涨落和探测器单元差异等原因，拟合得到的椭圆系列中心坐标有一定的偏差。通常会出现以下3种情况：①椭圆系中心坐标分布在很小的区域内，可以取这些点的重心作为最终的重建芯位坐标；②椭圆系中心坐标总体分布密集，个别偏离较远，需要去除偏离较远的那些坐标值，以密集点群的重心作为芯位(如图 5)；③椭圆系中心坐标分布极为离散，这种事例作为无效事例处理。decut是用来量化判断离散度的标准，通常decut的值越小，得到的芯位精度越高，但是太小的decut可能丢掉部分有效数据，导致重建效率很低。通过多次尝试比对，我们选择5 m这一数值，能够兼顾芯位精度与重建效率。在上述第2种情况中，通常中低能事例的小光电子数数据组拟合的椭圆中心坐标与芯位偏差较大，这是由于光电子数比较小的探测器单元信号信噪比低；而高能事例的大光电子数数据组拟合的椭圆中心坐标与芯位偏差较大，可能是因为高能事例芯位附近的光电倍增管已经饱和。
|
| 图 5 椭圆拟合法重建出的各个椭圆中心和真实芯位，其中蓝色+是各组数据重建得到的椭圆中心，红色*是真实芯位 Fig. 5 The center of each ellipse and the real core position reconstructed by the ellipse fitting method, where the blue "+" is the center of the ellipse reconstructed from each group of data, and the red "*" is the real core position |
(4) x0i和y0i的标准差。用x0i和y0i的标准差作为参数来判断这一系列椭圆中心的离散程度，若标准差很大，说明一些数据组几乎不包含芯位信息，这些光电子数主要是由噪声造成的，因此把这些数据当做噪声去除；而椭圆中心坐标很集中则表明与之对应的数据中噪声比很低，保留其对芯位坐标的贡献。图 5是一个事例分组重建的各个椭圆中心坐标，最终只保留红框内的密集点求平均作为最终的芯位重建结果。
(5) NBN与NBNmin。得到满足离散度要求的椭圆中心坐标平均值(x0, y0) 后，还需要利用附近10 m圆和10~20 m的环内光电子数总和之比NBN=N10/N20判断该点周围的光电子数分布符合芯位附近的光电子数分布特点，NBN > NBNmin表明重建芯位误差小于10 m，可以把计算结果作为重建芯位的最终结果，其中NBNmin的取值选择可以通过数值估计和数据统计得到。在需要重建阵列边界附近的芯外事例时，考虑到N10可能为0，可以用NBN=N20/N30或NBN=N30/N40也能达到同样的判断效果，只是涨落变大且精度变差，只能分别判断重建芯位误差是否低于20 m和30 m。
3 椭圆拟合法重建结果蒙特卡罗模拟是研究宇宙线物理的重要手段，本文借助Corsika模拟的广延大气簇射和Geant4模拟的探测器响应数据来验证椭圆拟合法重建芯位的实际效果。选取能量为10 TeV~10 PeV的质子事例，天顶角为0~40°，服从cos2θ分布，方位角为0~360°，服从均匀分布。探测器参考LHAASO-WCDA1，各项参数保持一致，在坐标系中探测器阵列区域为(-75 m < x < 75 m, -75 m < y < 75 m)，在数据中加入60 kHz的噪声，并利用遍历法取200 ns的时间窗口去噪之后信噪比约为70/22。我们分别用重心法和椭圆拟合法来重建事例并比较重建精度和事例利用率。首先需要验证本算法对于芯位位于探测器阵列边沿事例的重建能力，选取真实芯位(xct, yct)位于(75 m < x < 85 m, -75 m < y < 75 m) 探测器阵列外紧贴边界的能量为10~100 TeV质子事例重建芯位，得到结果如图 6。从图 6可以看出，重建芯位与真实芯位分布基本一致，从图 6 (b)可以看出, 大多数重建芯位的偏差
|
| 图 6 (a) 真实芯位和重建芯位分布；(b) 重建芯位偏差统计 Fig. 6 (a) Real core position and reconstructed core position distribution; (b) reconstructed core position deviation statistics |
为了定量描述该方法重建芯位的精度和对探测器面积的利用率，我们分别取4 000个能量为1~10 TeV, 10~100 TeV, 100~1 000 TeV和1~10 PeV，真实芯位均匀分布在(-150 m < x < 150 m, -150 m < y < 150 m) 区域的质子事例，用重心法和椭圆拟合法重建芯位。由于真实芯位均匀分布在300 m × 300 m的区域，而探测器阵列区域为150 m × 150 m，因此取重建成功的事例数与芯位位于探测器阵列的事例数之比作为探测器面积利用率
|
| 图 7 (a) 椭圆拟合法和重心法芯位重建分辨率；(b) 探测器阵列面积利用率 Fig. 7 (a) Ellipse fitting method and center of gravity method core position reconstruction resolution; (b) detector array area utilization |
从图 7可以看出，对于所有能量段的质子事例，椭圆拟合法芯位重建精度比重心法要好，并且椭圆拟合法能够显著提高探测器阵列利用率，随着能量增大，其增幅更加显著。对于1~10 TeV的事例，椭圆拟合法的探测器面积利用率略小于1，这意味着只有少数芯位落在探测器边沿的事例重建不好，面积利用率随着能量增大而增大，对于能量为1~10 PeV的事例，面积利用率大于1.3，这意味着约有30%位于探测器阵列外部事例的芯位重建符合精度要求。图 7 (a)椭圆拟合法重建芯位分辨率随着能量增大而增大的现象不符合预期，这是由于把大量芯外事例纳入统计量，导致面积利用率增大造成的，通过归一化面积利用率后的芯位分辨率可以得到证实。各能段面积利用率为0.67的情况下椭圆法和重心法的重建分辨率如图 8。在重建效率相同的情况下，椭圆拟合法的芯位分辨率显著优于重心法；并且两种方法的芯位分辨率都随着能量增加而降低，其中椭圆拟合法由于本身能够排除大部分噪声干扰，重建精度随能量变化不明显。
|
| 图 8 面积利用率为0.69时各能段的芯位重建分辨率 Fig. 8 Reconstruction resolution of core positions for each energy segment when the area utilization ratio is 0.69 |
通过蒙特卡罗模拟数据检验，椭圆拟合法相比重心法具有明显优势，椭圆拟合法的优势主要体现在两个方面：(1) 充分考虑广延大气簇射的对称性特征，利用对称性能够准确重建探测器阵列边沿乃至探测器区域之外的事例；(2) 通过分组重建的方式有效消除了噪声对芯位重建的影响，从而显著提高芯位重建的精度。从150 m × 150 m的WCDA-I模拟数据的重建结果来看，椭圆拟合法重建芯内事例，芯位重建精度可以达到5 m以下，各能段分辨率不到重心法的50%，重建精度相比于重心法有显著提高；另一方面能够解决重心法无法准确重建探测器阵列边沿事例芯位的问题，重建低能事例的面积利用率约比重心法高40%，重建高能事例的面积利用率则接近重心法的2倍；与重心法相比，其芯位重建精度和探测器阵列面积利用率都有显著改善。
椭圆拟合法的噪声抑制能力源自广延大气簇射次级粒子触发的NPE信号，等高线为一系列同心椭圆而噪声信号均匀分布，根据这一原理设计去噪算法，在遍历法去噪的基础上利用芯位信息实现深度去噪。此外，前锋面次级粒子到达时间关于簇芯对称分布，我们可以基于这一特征通过分区重建方向规避锥面修正参数依赖带来的问题，并且利用方向信息实现深度去噪。
| [1] | CAO Z, CHEN M J, CHEN S Z, et al. Introduction to Large High Altitude Air Shower Observatory (LHAASO)[J]. Chinese Astronomy and Astrophysics, 2019, 43(4): 457–478. DOI: 10.1016/j.chinastron.2019.11.001 |
| [2] |
祝凤荣. 用ARGO-YBJ实验的SPT寻找~GeV的伽玛暴[M]. 成都: 西南交通大学出版社, 2013. ZHU F R. Searching for gamma-ray bursts of~GeV with SPT of ARGO-YBJ experiment[M]. Chengdu: Southwest Jiaotong University Press, 2013. |
| [3] |
陈宝民, 高博, 陈明君, 等. WCDA时间标定系统的研制与测试[J]. 天文研究与技术, 2020, 17(3): 399–404 CHEN B M, GAO B, CHEN M J, et al. Development and testing of WCDA time calibration system[J]. Astronomical Research & Technology, 2020, 17(3): 399–404. DOI: 10.14005/j.cnki.issn1672-7673.20191129.001 |
| [4] | JOSHI V, HINTON J, SCHOORLEMMER H, et al. A template-based γ-ray reconstruction method for air shower arrays[J]. Journal of Cosmology and Astroparticle Physics, 2019, 2019(1): 012. DOI: 10.1088/1475-7516/2019/01/012 |
| [5] | FENG Y L, ZHANG Y, CHEN T L, et al. Lateral distribution of EAS muons measured for the primary cosmic ray energy around 100 TeV[J]. Chinese Physics C, 2019, 43(7): 97–102. |
| [6] | GREISEN K. Progress in cosmic ray physics VOL. 3[M]. [S: l]: North-Holland Publishing, 1956. |
| [7] |
吴文雄, 左雄, 肖刚, 等. LHAASO-MD光电倍增管的性能测试[J]. 天文研究与技术, 2020, 17(2): 258–264 WU W X, ZUO X, XIAO G, et al. The performance test of LHAASO-MD photomultiplier[J]. Astronomical Research & Technology, 2020, 17(2): 258–264. DOI: 10.14005/j.cnki.issn1672-7673.20191112.005 |
| [8] | ABU-ZAYYAD T, BELOV K, BIRD D, et al. A measurement of the average longitudinal development profile of cosmic ray air showers between 1017 and 1018 eV[J]. Astroparticle Physics, 2001, 16: 1–11. DOI: 10.1016/S0927-6505(00)00170-5 |
| [9] | COLLABORATION H, ABEYSEKARA A U, ALFARO R, et al. The HAWC Gamma-Ray Observatory: design, calibration, and operation[C]//Proceedings of the 33rd International Cosmic Ray Conference. 2013. |
| [10] | ANTONOV R A, BONVECH E A, CHERNOV D V, et al. Spatial and temporal structure of EAS reflected Cherenkov light signal[J]. Astroparticle Physics, 2019, 108: 24–39. DOI: 10.1016/j.astropartphys.2019.01.002 |
| [11] | TOMA G. Investigation of the EAS lateral particle density at 500 m distance from shower core[C]//Proceedings of the AIP Conference. 2008: 549-553. |
| [12] | BARTOLI B. EAS age determination from the study of the lateral distribution of charged particles near the shower axis with the ARGO-YBJ experiment[J]. Astroparticle Physics, 2017, 93: 46–55. DOI: 10.1016/j.astropartphys.2017.06.003 |
| [13] | LI C, HE H H, XIAO G, et al. Measurement of muonic and electromagnetic components in cosmic ray air showers using LHAASO-KM2A prototype array[J]. Physical Review D, 2018, 98(4): 042001. DOI: 10.1103/PhysRevD.98.042001 |
| [14] | BRÂNCUS I M. Correlated features of arrival time and angle-of-incidence distributions of EAS muons[J]. Astroparticle Physics, 1997, 7(4): 343–356. DOI: 10.1016/S0927-6505(97)00030-3 |
| [15] | 王晓洁. LHAASO-WCDA本底噪声过滤与簇射事例重建方法的研究[D]. 北京: 中国科学院大学, 2017. WANG X J. Research on background noise filtering and shower case reconstruction method of LHAASO-WCDA[D]. Beijing: Unversity of Chinese Academy of Sciences, 2017. |
| [16] | HEDAYATI H, MORADI A, EMAMI M. A statistical method for reconstructing the core location of an extensive air shower[J]. The Astrophysical Journal, 2015, 810(1): 68. DOI: 10.1088/0004-637X/810/1/68 |
| [17] | CHEN S Z. Status of LHAASO official software: simulation and reconstruction[C]//Proceedings of the 8th International Workshop on Air Shower Detection At High Altitudes. 2017: 1-37. |
| [18] | RAVIGNANI D, SUPANITSKY A D, MELO D. Reconstruction of air shower muon densities using segmented counters with time resolution[J]. Astroparticle Physics, 2016, 82: 108–116. DOI: 10.1016/j.astropartphys.2016.06.001 |
| [19] | APEL W D, ARTEAGA-VELAZQUEZ J C, BEKK K, et al. Cosmic ray energy reconstruction from the S (500) observable recorded in the KASCADE-Grande air shower experiment[J]. Astroparticle Physics, 2016, 77: 21–31. DOI: 10.1016/j.astropartphys.2015.12.002 |
| [20] | CAO Z. A future project at tibet: the large high altitude air shower observatory (LHAASO)[J]. Chinese Physics C, 2010, 34(2): 249–252. DOI: 10.1088/1674-1137/34/2/018 |
| [21] |
曹臻, 陈明君, 陈松战, 等. 高海拔宇宙线观测站LHAASO概况[J]. 天文学报, 2019, 60(3): 3–18 CAO Z, CHEN M J, CHEN S Z, et al. Introduction to Large High Altitude Air Shower Observatory (LHAASO)[J]. Acta Astronomica Sinica, 2019, 60(3): 3–18. |
| [22] | DEY R K, DAM S, BASAK A. A novel approach for deducing the mass composition of cosmic rays from lateral densities of EAS particles[J]. Europhysics Letters, 2019, 127(3): 39002. DOI: 10.1209/0295-5075/127/39002 |