# Automated Chest X-Ray Report Generator Using Multi-Model Deep Learning
Approach

1st Arief Purnama Muharram  School of Electrical Engineering and Informatics
Institut Teknologi Bandung
Bandung, Indonesia
<EMAIL_ADDRESS>2nd Hollyana Puteri Haryono  School of Electrical
Engineering and Informatics
Institut Teknologi Bandung
Bandung, Indonesia
<EMAIL_ADDRESS>3rd Abassi Haji Juma  School of Electrical
Engineering and Informatics
Institut Teknologi Bandung
Bandung, Indonesia
<EMAIL_ADDRESS>4th Ira Puspasari  Computer Engineering
Universitas Dinamika
Surabaya, Indonesia
<EMAIL_ADDRESS>5th Nugraha Priya Utama* *Corresponding author School of
Electrical Engineering and Informatics
Institut Teknologi Bandung
Bandung, Indonesia
<EMAIL_ADDRESS>
###### Abstract

Reading and interpreting chest X-ray images is one of the most radiologist’s
routines. However, it still can be challenging, even for the most experienced
ones. Therefore, we proposed a multi-model deep learning-based automated chest
X-ray report generator system designed to assist radiologists in their work.
The basic idea of the proposed system is by utilizing multi binary-
classification models for detecting multi abnormalities, with each model
responsible for detecting one abnormality, in a single image. In this study,
we limited the radiology abnormalities detection to only cardiomegaly, lung
effusion, and consolidation. The system generates a radiology report by
performing the following three steps: image pre-processing, utilizing deep
learning models to detect abnormalities, and producing a report. The aim of
the image pre-processing step is to standardize the input by scaling it to
128x128 pixels and slicing it into three segments, which covers the upper,
lower, and middle parts of the lung. After pre-processing, each corresponding
model classifies the image, resulting in a 0 (zero) for no abnormality
detected and a 1 (one) for the presence of an abnormality. The prediction
outputs of each model are then concatenated to form a ’result code’. The
’result code’ is used to construct a report by selecting the appropriate pre-
determined sentence for each detected abnormality in the report generation
step. The proposed system is expected to reduce the workload of radiologists
and increase the accuracy of chest X-ray diagnosis.

###### Index Terms:

chest x-ray, radiology, medical report, multi-model, deep learning

##  I Introduction

Deep learning for image processing has rapidly advanced in recent years [1],
demonstrating promising results across a variety of domains, including
radiology [2]. These kind of advancements have enabled improvements in
radiological diagnosis and patient care[2]. Therefore, there is a demand for
automated systems that can accurately and quickly process the growing volume
of medical images [2]. Chest X-ray (CXR) imaging, in particular, is one of the
most often utilized radiologic diagnostic tools, although its interpretation
can be difficult even for the experienced radiologists[3]. Radiologists can
benefit from automated CXR report generation since it allows them to focus on
more difficult cases while reducing the possibility of errors.

In this study, we proposed a multi-model deep learning-based automated CXR
report generator system to help radiologists’ work better and faster. The
proposed system is designed to detect abnormalities in a CXR image and then
produce the corresponding report. Our approach was motivated by previous
research that successfully applied deep learning algorithms to the processing
of chest radiographs [4, 5, 6, 7, 8, 9, 10]. However, most of the studies
mentioned have primarily focused on detecting or classifying a single
abnormality in the image. Therefore, in this study, we aimed to detect
multiple abnormalities in radiology images. We hypothesized that using multi-
model, with each model responsible for an abnormal finding, can detect
multiple abnormalities present within a single image.

This paper will be structured as follows: Firstly, we will discuss related
work on deep learning in radiology and mention several previous research
studies in the field. Next, we will present the methodology and technical
aspects of our suggested system. In the Results section, we will describe and
discuss our experimental findings. Finally, we will conclude our study.

TABLE I: Literature Review Author | Year | Aim | Algorithm Used |  Performance
---|---|---|---|---
Rajkomar et al. [4] | 2017 |  Developing a deep convolutional neural network for classifying the view orientations of chest radiographs. |  GoogLeNet |  Accuracy of 100% and 38 image classifications per second.
Lakhani et al. [6] | 2017 |  Developing CNNs for automated classification of pulmonary tuberculosis using chest radiograph |  AlexNet, GoogLeNet |  The best performance is achieved by the ensemble of AlexNet and GoogLeNet, with an Area Under the Curve (AUC) of 0.99.
Cicero et al. [5] | 2017 |  Developing and test a deep convolutional neural network for frontal chest radiograph abnormality classification and computer-aided detection |  GoogLeNet |  The sensitivity, specificity, and AUC for pleural effusion were 0.91, 0.91, and 0.962, respectively; for pulmonary edema, they were 0.82, 0.82, and 0.868; for consolidation, 0.74, 0.75, and 0.850; for cardiomegaly, 0.81, 0.80, and 0.875; and for pneumothorax, 0.78, 0.78, and 0.861.
Saha et al. [7] | 2020 |  Developing multi model based ensemble to detect COVID-19 from CXR Images |  Multiple CNN models (CNNs, MobileNet, InceptionV3, DenseNet201, DenseNet121 and Xception) ensembled using feature concatenation and decision fusion |  Accuracy of 96% using decision fusion scheme for 3 class classification and 89.21% for 4 class classification. Accuracy of 95.84% using feature concatenation scheme for 3 class classification and 89.26% for 4 class classification
Deb et al. [8] | 2022 |  Developing a multi model ensemble based DCNN structure |  VGGNet, GoogLeNet, DenseNet, and NASNet |  An accuracy of 88.98% was achieved for the three-class classification, and a 98.58% accuracy was attained for the binary class classification using the publicly available dataset.
Tang et al. [9] | 2020 |  Developing deep CNN model to classify normal vs abnormal CXRs |  AlexNet, VGG, GoogLeNet, ResNet, and DenseNet |  ROC score for ResNet18, ResNet50, GoogLeNet, and DenseNet are by around 98% using both NIH CXR and WCMC pediatric dataset.
Bhatt et al. [10] | 2021 |  Developing classification model for pulmonary consolidations detection in CXR images using transfer learning and progressive resizing techniques |  VGG-19 (baseline for binary classification) and EfficientNet-B3 (baseline for multi-class classification) |  Achieved 100% accuracy
Amjoud et al. [11] | 2021 |  Developing automated report generation for CXR using Transformer-based model |  Transfer learning and transformer-based approach |  Achieved BLUE-1, BLUE-2, and ROUGE metrics, where the scores were 0.479, 0.359, and 0.380, respectively.
Ghadekar et al. [12] | 2021 |  Developing report generation using VGG for chest diseases diagnosis and encoder-decoder model for report generation. |  VGG16 and encoder-decoder model |  Accuracy of 88% for classification images and 96% for report generation.

##  II Related work

In recent years, medical image analysis has emerged as a major use of deep
learning, with the potential to aid radiologists in diagnosing a disease. Deep
learning algorithms, in particular, have shown promise in the identification
and classification of abnormalities on chest radiographs. One of the popular
and widely-used deep learning algorithms in the field is Convolutional Neural
Network (CNN).

LeCun et al. [1] present an in-depth overview of deep learning, including
CNNs, which have been demonstrated to be successful for image classification
tasks. CNN is a form of artificial neural network (ANN) that is extensively
employed in image recognition and processing applications. It is based on the
structure and function of the visual cortex in animals, which is responsible
for visual information processing [1].

In a CNN, input data is fed into a sequence of convolutional layers, which
extract features from the input image by applying a set of filters. To inject
non-linearity into the model, the output of each convolutional layer is then
processed through a non-linear activation function such as ReLU (Rectified
Linear Unit). As a result, the input data is represented by a set of high-
level features that become increasingly complex and abstract.

The output of the convolutional layers is often sent into one or more fully
connected layers, which conduct classification or regression based on the
retrieved features. The network’s ultimate output is a probability
distribution over the various classes or values of the output variable.

CNNs have proven to be extremely effective in image classification, object
detection, and segmentation tasks, and have been utilized in a variety of
applications including as self-driving cars, medical diagnostics, and facial
recognition.

CNNs were shown by Rajkomar et al. [4] to be capable of classifying view
orientations of chest radiographs with excellent accuracy. Lakhani et al. [6]
used chest radiographs to construct a CNN-based model for the automated
classification of pulmonary tuberculosis, obtaining high performance and
indicating the promise for deep learning in the disease detection. Cicero et
al. [5] used a CNN to detect and classify abnormalities on chest radiographs,
with good sensitivity and specificity. Tang et al. [9] compare several
existing CNN architectures such as AlexNet, VGGNet, ResNet, GoogLeNet, and
DenseNet with tested on NIH CXR dataset using ROC score for evaluation. It
found that DenseNet, GoogLeNet, and ResNet has performance of ROC score more
than 98%. Bhatt et al. [10] use modified VGG-19 and EfficientNet-B3
architecture and applied transfer learning technique to identify pneumonia and
COVID-19 in CXR images. Table I describes in detail previous research on deep
learning in radiology.

There are several papers that used multi-model based ensembling approaches for
detection of COVID19 from CXR images. Saha et al. [7] use ensembling
approaches with both feature concatenation and decision fusion technique. Deb
et al. [8] develop a multi-model ensemble based on DCNN architecture with
combining VGGnet, GoogLeNet, DenseNet, and NASNet.

Report Generation of CXR using deep learning approach has been done before by
Amjoud et al. [11] and Ghadekar et al. [12]. Encoder-Decoder approach has been
the focus in Ghadekar et al. [12] and giving a results of 96% of accuracy.
However, they do not mentioned the encoder-decoder architecture used, nor do
they explain the characteristics of the generated report. Amjoud et al. [11]
focused on combining CNN for CXR image classification and transformers for
report generation.

Compared to the mentioned studies, our proposed method utilizes a multi-model
approach to detect multiple abnormalities present in the image. The results
from each model are then used to generate the report. We argue that using
multiple models, each dedicated to a specific abnormality, leads to more
’responsible’ results.

##  III Methodology

###  III-A System Design

Figure 1: Three fundamental steps in our proposed methodology

Based on our literature review, we are motivated to create a deep-learning
model for quantifying and localizing abnormalities on chest radiographs and
generating corresponding X-ray reports, to support radiologists in diagnosing
a disease. Our proposed methodology consists of three fundamental steps: image
pre-processing, CXR deep learning models, and report production (Figure 1).

Image pre-processing is the first step of the entire process, as illustrated
in Figure 2. In this step, each image is resized to 128x128 and square-cropped
to ensure consistent dimensionality. Then, the image is divided into three
segments (Segment I, II, and III) to reduce the amount of ”information” passed
to the deep learning model in the subsequent step. Segment I covers the upper
part of the lung, segment II covers the lower part of the lung, and segment
III covers the middle part of the lung.

Figure 2: Image pre-processing step

The second step involves training deep learning models to detect each
abnormality. Three models were trained, one for detecting cardiomegaly,
another for detecting lung effusion, and the third for detecting
consolidation. However, instead of using the whole image as input for the
model, we used the corresponding sliced image based on the specific
abnormality of interest (Figure 3). For example, we used segment II as the
input for the cardiomegaly and lung effusion models because these
abnormalities are more likely to be found in the lower part of the lungs,
while segment III was used as input for the consolidation model because this
abnormality is more likely to be found in the middle part of the lungs. The
NIH CXR dataset was used to train the models [15]. For the experiments, we
used the ”one-factor-at-a-time” strategy, using the learning rate, optimizer,
and pre-trained models as the observable factors. In addition, we compared the
ResNet18 [13], ResNet50 [13], and GoogleLeNet [14] pre-trained models, which
are popular CNN algorithms.

Figure 3: Multi-models breakdown and result code

We trained our models on a high-performance computing server with an Intel
Xeon Silver 4208 CPU 3.2 GHz, 256 GB of RAM, and an NVIDIA Quadro RTX 5000 GPU
with 16 GB of RAM. For our deep learning library, we used Torch version 1.13.1
and Torchvision version 0.14.1 in a Python 3.9.6 environment, and for
deployment, we used Gradio version 3.28.3.

Finally, the last step is report generation. The purpose of this step is to
map the result code, which is an aggregation of each model’s prediction
output, to the corresponding sentence in the master text (Figure 4). For
example, if a cardiomegaly abnormality is detected (coded as “1”), the
sentence “Terdapat kardiomegali, CTR $<$ 50%” will be selected from the master
text. Otherwise, the sentence “Bentuk jantung baik, tidak ditemukan
kardiomegali” will be used.

Figure 4:  Report generation is based on the result code and master text

###  III-B Evaluation Strategy

We used accuracy as our evaluation metric, both for individual trained models
and for overall system performance. Accuracy is determined by calculating the
total number of correct predictions divided by the total number of
predictions. However, when assessing the system’s performance, a true
prediction is only achieved when all constituent models provide accurate
predictions. For instance, if the prediction labels for the cardiomegaly, lung
effusion, and consolidation models are [0, 0, 1], respectively, and the ground
truth labels are [0, 0, 0], then the prediction result is considered
incorrect.

###  III-C Problem Limitations

There are several limitations in our study that should be acknowledged,
including:

  * •

Radiology image limitation: The suggested method was trained using adult
frontal CXRs and is only intended to be used for adult patients’ frontal CXR
images. Therefore, the system may not be suitable for interpreting X-rays of
pediatric patients.

  * •

Abnormalities detection limitation: The method is designed to identify
specific abnormalities, namely cardiomegaly, lung effusion, and consolidation.
As a result, any additional abnormalities present in the image that fall
outside the scope of this study are not be detected by the system.

##  IV Result

###  IV-A Dataset Pre-processing

We used The NIH CXR Dataset [15] to train our models, which comprises hundreds
of thousands of CXR images with 14 different labels for abnormality findings.
However, we only used 3 of the 14 abnormality labels, namely cardiomegaly,
lung effusion, and consolidation. To ensure data quality, we employed domain
experts to conduct reannotation. For each abnormality, we selected 2,000
images at random (1,000 images with the presence of abnormality and 1,000
without). Any images that did not meet our criteria, such as low image
quality, incorrect orientation, or misdiagnosis, were excluded. Following
annotation, the dataset was randomly split into the training (70%) and testing
(30%) datasets. Table II presents the results of our annotated dataset.

TABLE II: Total number of images in the reannotated dataset Abnormalities | Labels | Total (N)
---|---|---
False (n) | True (n)
Cardiomegaly | 615 | 793 | 1408
Lung Effusion | 658 | 896 | 1554
Consolidation | 454 | 652 | 1106

To evaluate the system’s performance, we created a separate dataset, named the
”system evaluation dataset,” by randomly selecting 200 images from The NIH CXR
dataset. In contrast to the model training dataset, this dataset comprises
three labels: cardiomegaly, lung effusion, and consolidation, each with a
binary value of 0 for absence or 1 for presence of abnormality. This dataset
is designed to assess the system’s ability to accurately predict all three
labels simultaneously. Additionally, domain experts have reannotated the
dataset to ensure its quality, just as we did with the model training dataset.

###  IV-B Modeling

There were several parameters that we considered for hyperparameter tuning,
such as learning rate, optimizer, and pretrained model. Our hyperparameter
tuning strategy was based on the ”one-factor-at-a-time” approach, where we
varied one parameter at a time while keeping the others at a default value.
For each parameter value, we compared the resulting training and testing
accuracy. Hyperparameter tuning was performed for each model, and the results
are presented in Tables III, IV, and V.

After completing hyperparameter tuning, we found that the optimal values were
different for each abnormality. The most optimal learning rates were 1e-3 for
effusion and cardiomegaly, and 1e-4 for consolidation. Adam was the best
optimizer since we found that SGD did not produce the best accuracy. Also, the
most optimal pre-trained model for each abnormalities were different, such as
GoogLeNet for cardiomegaly, ResNet50 for effusion, and ResNet18 for
consolidation. Using these optimal parameters for each abnormality, we trained
the final model and recorded the training and testing accuracy as shown on
Table VI. The final models will be the models used in the system.

TABLE III: Hyperparameter Tuning Result for Cardiomegaly Step | Hyperparameter | Value | Accuracy
---|---|---|---
Training | Testing
1 | Learning Rate | 1e-2 | 75.4% | 80.5%
1e-3 (default) | 87.1% | 87%
1e-4 | 82.6% | 81.5%
1e-5 | 78.1% | 79.8%
2 | Optimizer |  Adam (default) | 87.1% | 87%
SGD | 72.4% | 80.5%
3 | Pre-trained Model | ResNet50 (default) | 87.1% | 87%
ResNet18 | 85% | 90.3%
GoogLeNet | 90.1% | 90.1%
TABLE IV: Hyperparameter Tuning Result for Lung Effusion Step | Hyperparameter | Value | Accuracy
---|---|---|---
Training | Testing
1 | Learning Rate | 1e-2 | 77.5% | 78.4%
1e-3 (default) | 90.2% | 91.8%
1e-4 | 87% | 87%
1e-5 | 82% | 84.6%
2 | Optimizer |  Adam (default) | 90.2% | 91.8%
SGD | 76.4% | 79.4%
3 | Pre-trained Model |  ResNet50 (default) | 90.2% | 91.8%
ResNet18 | 91.5% | 87.4%
GoogLeNet | 87.8% | 89.4%
TABLE V: Hyperparameter Tuning Result for Consolidation Step | Hyperparameter | Value | Accuracy
---|---|---|---
Training | Testing
1 | Learning Rate | 1e-2 | 67.5% | 68.9%
1e-3 (default) | 79.1% | 78%
1e-4 | 79.1% | 81%
1e-5 | 69.6% | 74.7%
2 | Optimizer |  Adam (default) | 79.1% | 78%
SGD | 66% | 65.6%
3 | Pre-trained Model | ResNet50 (default) | 79.1% | 78%
ResNet18 | 82.7% | 78.4%
GoogLeNet | 78.5% | 73.3%
TABLE VI: Result of each final model Abnormality | Accuracy
---|---
Training | Testing
Cardiomegaly | 89.2% | 88.7%
Lung Effusion | 88.6% | 86.1%
Consolidation | 75.9% | 81%

###  IV-C System Evaluation

Figure 5: System example TABLE VII: System Evaluation Abnormality | Accuracy
---|---
Cardiomegaly | 52%
Lung Effusion | 80%
Consolidation | 40%
Cardiomegaly + Lung Effusion + Consolidation | 20%

We successfully created and deployed the system in the form of a web-based
application (Figure 5). Tthe system’s performance evaluation was done using
the system evaluation dataset. For a correct prediction, the system needed to
accurately predict all three abnormalities; any incorrect prediction would
render the overall result incorrect. Despite our hypothesis that the multi-
model approach would yield promising results, the system’s accuracy was only
20% (Table VII).

Our proposed system exhibits significantly lower performance in comparison to
the studies conducted by Amjoud et al. [11] and Ghadekar et al. [12] in the
generation of radiology reports. To investigate the poor performance of the
system, we conducted an error analysis that included a model breakdown
analysis. We discovered that the accuracy of the cardiomegaly, lung effusion,
and consolidation models towards the system evaluation dataset were only 52%,
80%, and 40%, respectively (Table VII). Of the three models, the worst
performance results were obtained by the cardiomegaly and consolidation
models.

The poor results obtained by the cardiomegaly and consolidation models
indicate that these models failed to generalize well. Although the system
evaluation dataset and model training dataset were sourced from the same
dataset (the NIH CXR Dataset [15]), this doesn’t guarantee that the data is
the ”same”. Medical images are highly heterogeneous, which can increase the
risk of overfitting and lead to a loss of generalizability for the models [2].

Furthermore, the errors of each model contribute to the poor performance of
the system. Based on the accuracy, it is evident that the error rates of the
cardiomegaly, lung effusion, and consolidation models were 48% (1), 20% (2),
and 60% (3), respectively. The combination of these errors could make it
difficult for the system to make accurate predictions, and this could explain
the poor performance results.

Below is our mathematical proof of probability to explain the poor performance
of our system. Let $A$, $B$, and $C$ denote the cardiomegaly, lung effusion,
and consolidation models, respectively. $P_{correct}(A)$, $P_{correct}(B)$,
and $P_{correct}(C)$ represent the probabilities of accurately predicting the
cardiomegaly, lung effusion, and consolidation models, respectively. In this
case, the models’ probabilities of accurate prediction are the models’
accuracies themselves. Since the models function independently, we can assume
that the events are also independent. We can calculate the probabilities of
correct and error predictions using the following formula.

$P_{error}(A)=1-P_{correct}(A)=1-0.52=0.48$ (1)

$P_{error}(B)=1-P_{correct}(B)=1-0.8=0.2$ (2)

$P_{error}(C)=1-P_{correct}(C)=1-0.4=0.6$ (3)

$P_{correct}(A\cap B\cap C)=P_{correct}(A)\times P_{correct}(B)\times
P_{correct}(C)=0.52\times 0.8\times 0.4=0.1664$ (4)

$P_{error}(A\cup B\cup
C)=P_{error}(A)+P_{error}(B)+P_{error}(C)-P_{error}(A\cap B)-P_{error}(A\cap
C)-P_{error}(B\cap C)+P_{error}(A\cap B\cap C)=0.48+0.2+0.6-(0.48\times
0.2)-(0.48\times 0.6)-(0.2\times 0.6)+(0.48\times 0.2\times 0.6)=0.8336$ (5)

The calculation shows that the probability of accurate predictions can be as
low as $\approx 16\%$ while the probability of error predictions can be as
high as $\approx 83\%$, which explains the system’s poor performance in
mathematically. Additionally, this highlights that the performance of the
overall system is affected by the individual models’ performance, which could
explain the system’s inferior performance when compared to the Amjoud et al.
[11] and Ghadekar et al. [12] studies.

##  V Conclusion

We developed an automated CXR report generator based on multi-model deep
learning, capable of detecting and classifying abnormalities in CXR images.
Unfortunately, our system’s overall accuracy was only 20%, indicating that our
multi-model technique of using separate models for each abnormality does not
have the potential to significantly improve the accuracy and efficiency of
radiological report generation.

Nonetheless, our method has a number of limitations that require further
investigation in future research. One of the limitations is the dataset used
in this study, which is unable to capture the heterogeneity of medical images,
thereby limiting the model’s ability to generalize the data. Consequently, the
model’s underperformance resulted in poor system performance. Future studies
should aim to increase the amount of validated radiology data or develop
alternative methods or algorithms to detect multiple abnormalities in a single
radiology image.

## Acknowledgment

This work was the IF5200 Applied Research Project course project of the first,
second, and third authors, under the supervision of the fourth and fifth
authors. We would like to express our gratitude to the Informatics Master
Study Program, School of Electrical Engineering and Informatics, Institut
Teknologi Bandung, for providing us access to their high-performance computing
server. Our research project would not have been possible without this vital
resource. Furthermore, we also extend our appreciation to our colleague,
Komang Shary Karismaputri, M.D., a radiology resident from Cipto Mangunkusumo
Hospital at the Faculty of Medicine, Universitas Indonesia. Her insightful
contributions to this research, particularly in providing valuable insights
into the interpretation of radiology images, have greatly enriched the quality
and depth of our study.

## References

  * [1] Y. LeCun, Y. Bengio, and G. Hinton, ”Deep learning,” Nature, vol. 521, pp. 436-444, 2015.
  * [2] J. H. Thrall, X. Li, Q. Li, C. Cruz, S. Do, K. Dreyer, and J. Brink, ”Artificial Intelligence and Machine Learning in Radiology: Opportunities, Challenges, Pitfalls, and Criteria for Success,” J. Am. Coll. Radiol., vol. 15, no. 3, pp. 504-508, Mar. 2018.
  * [3] C. S. Lee, P. G. Nagy, S. J. Weaver, and D. E. Newman-Toker, ”Cognitive and system factors contributing to diagnostic errors in radiology,” AJR. American journal of roentgenology, vol. 201, no. 3, pp. 611-617, 2013.
  * [4] A. Rajkomar, S. Lingam, A. G. Taylor, M. Blum, and J. Mongan, ”High-throughput classification of radiographs using deep convolutional neural networks,” J. Digital Imaging, vol. 30, no. 1, pp. 95-101, Feb. 2017.
  * [5] M. Cicero, A. Bilbily, E. Colak, D. Kontos, and J. Mermelstein, ”Training and validating a deep convolutional neural network for computer-aided detection and classification of abnormalities on frontal chest radiographs,” Investigative Radiology, vol. 52, no. 5, pp. 281-287, May 2017.
  * [6] P. Lakhani and B. Sundaram, ”Deep learning at chest radiography: automated classification of pulmonary tuberculosis by using convolutional neural networks,” Radiology, vol. 284, no. 2, pp. 574-582, Aug. 2017.
  * [7] O. Saha, J. Tasnim, T. Raihan, T. Mahmud, I. Ahmmed, and S. A. Fattah, ”A Multi-Model Based Ensembling Approach to Detect COVID-19 from Chest X-Ray Images,” in 2020 IEEE Region 10 Conference (TENCON).
  * [8] S. D. Deb, R. K. Jha, K. Jha, and P. S. Tripathi, ”A multi-model ensemble based deep convolution neural network structure for detection of COVID19,” Biomed. Signal Process. Control, vol. 71, p. 103126, Mar. 2022.
  * [9] Y. X. Tang, Y. B. Tang, Y. Peng, et al., ”Automated abnormality classification of chest radiographs using deep convolutional neural networks,” npj Digit. Med., vol. 3, p. 70, May 2020.
  * [10] A. Bhatt, A. Ganatra, and K. Kotecha, ”COVID-19 pulmonary consolidations detection in chest X-ray using progressive resizing and transfer learning techniques,” Heliyon, vol. 7, no. 6, p. e07211, Jun. 2021.
  * [11] A. B. Amjoud and M. Amrouch, ”Automatic Generation of Chest X-ray Reports Using a Transformer-based Deep Learning Model,” in 2021 Fifth International Conference On Intelligent Computing in Data Sciences (ICDS).
  * [12] P. Ghadekar, A. Jevrani, S. Dumpala, S. Dass, A. Pandey, and R. Bansal, ”Deep Learning model for diagnosis and report generation of lethal chest diseases using X-rays,” in 2021 International Conference on Artificial Intelligence and Machine Vision (AIMV), pp. 1-5.
  * [13] K. He, X. Zhang, S. Ren, and J. Sun, ”Deep Residual Learning for Image Recognition,” in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 770-778.
  * [14] C. Szegedy et al., ”Going deeper with convolutions,” in 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 1-9.
  * [15] X. Wang, Y. Peng, L. Lu, Z. Lu, M. Bagheri, and R. M. Summers, ”NIH Chest X-rays [Data set],” [Online]. Available: https://nihcc.app.box.com/v/ChestXray-NIHCC.


# A Double Layered Water Cherenkov Detector Array for Gamma-Ray Astronomy

Samridha Kunwar  Hazal Goksu  Jim Hinton  Harm Schoorlemmer  Andrew Smith
Werner Hofmann  Felix Werner

###### Abstract

Ground-level particle detection is now a well-established approach to TeV
$\gamma$-ray astronomy. Detection of Cherenkov light produced in water-filled
detection units is a proven and cost-effective approach. Here we discuss the
optimization of the units towards the future Southern Wide-field Gamma-ray
Observatory (SWGO). We investigate a configuration in which each Water
Cherenkov Detector (WCD) unit in the array comprises two chambers with black
or reflective walls and a single PMT in each chamber. A shallow lower chamber
with a PMT facing downwards enables muon tagging and the identification of
hadron-induced air showers, which are the primary source of background in
$\gamma$-ray astronomy. We investigate how $\gamma$/hadron separation power
and achievable angular resolution depend on the geometry and wall reflectivity
of the detector units in this configuration. We find that excellent angular
resolution, background rejection power and low-energy response are achievable
in this double-layer configuration, with some reflective surfaces in both
chambers.

###### keywords:

Gamma-Ray , Water Cherenkov Detector (WCD) Array , Simulations

††journal: Nucl. Instr. Meth. Phys. Res.

[inst1]organization=Max-Planck-Institut fur
Kernphysik,addressline=Saupfercheckweg 1, city=Heidelber, postcode=69117,
country=Germany

[inst3]organization=Department of Physics, University of Maryland,city=College
Park, postcode=20742, state=Maryland, country=USA

[inst2]organization=IMAPP, Radboud University Nijmegen,city=Nijmegen,
country=The Netherlands

##  1 Introduction

Ground-level particle-based detection of air showers is a rapidly developing
approach to $\gamma$-ray astronomy at very high energies, with cosmic-ray
protons and nuclei as the main source of background. High-density arrays
maximise the number of detectable particles per air shower and water-Cherenkov
based detector units have proven to be an effective way of achieving large
array area, as demonstrated by HAWC (High-Altitude Water Cherenkov) Abeysekara
et al. [1] and LHAASO (Large High Altitude Air Shower Observatory) coverage
Aharonian et al. [2]. The Southern Wide-field Gamma-ray Observatory (SWGO)
Abreu et al. [3], Hinton [4], is a project towards constructing a large
detector array in the southern hemisphere with an advanced detector design and
superior performance compared to both HAWC and LHAASO. This paper investigates
the reference design; a double layered WCD. Several other advanced designs
such as a shallow WCD with 4 PMT’s are being considered Conceição, R. et al.
[5]. The performance of WCD arrays is largely driven by high altitude Sinnis
[6], large array area and large fill-factor Schoorlemmer et al. [7], but the
particle detection thresholds of the individual WCD units and their response
characteristics will influence the threshold energy and performance of the
array (see e.g. Assis et al. [8]).

The role of an individual detector unit in a ground-particle-based
$\gamma$-ray instrument is to measure the local shower particle number or
energy density, assign a local arrival time and ideally provide additional
information for $\gamma$/hadron separation. Muon identification, as a means of
hadronic background rejection (see, e.g. Schoorlemmer et al. [7]), can be
implemented either using separate detector elements, such as the LHAASO buried
muon detectors Aharonian et al. [2] or by discrimination within a standard WCD
unit. Additionally, the topology of the shower amplitude distribution such as
the charge distribution close to the shower core and the clustering
characteristics of hits far from the core region also provide information for
effective $\gamma$/hadron discrimination as demonstrated by Milagro Atkins and
et al [9].

The two-chamber concept, comprising a water volume separated into optically
isolated top and bottom chambers, explored here, provides a (potentially)
cost-effective general-purpose element with reasonable time and amplitude
resolution as well as muon identification Letessier-Selvon et al. [10]. The
principle is well established: most shower muons will pass through the entire
detector element with only ionisation losses and produce Cherenkov light in
both chambers. Electrons and $\gamma$-rays generate cascades that penetrate to
the lower chamber only in the case of high energy initial particles, and in
any case, the comparison of the signals in both sections allows discrimination
from a through-going muon. As the lower chamber captures the cascade products
of only high energy shower particles, it can also extend the detector dynamic
range close to the shower core in the case of high energy primary
$\gamma$-rays.

The quantities to be optimised are the dimensions of the two chambers, the
reflective properties of the internal surfaces and the photosensors.

We first discuss particle detection efficiency, photon timing distributions
and background trigger rates of individual WCD units. The main current
instruments of this type, HAWC and the LHAASO WCD array Aharonian et al. [2],
employ non-reflective surfaces for detector water volumes. This approach
reduces the time spread of light
reaching the photo-sensor but results in a less-uniform response and reduced
overall detection efficiency compared to reflective surfaces.

After a brief overview of the double layered candidate for the SWGO WCD design
(Section 2) and of the simulation tools used (Section 3), we explore the
trade-offs associated with the reflectivity of the unit surfaces and the
overall unit dimensions (Sections 4 and 5). In these sections we also compare
our results from the double layered WCDs with single-chamber WCDs similar to
HAWC and LHAASO units. In a second step we then address the energy threshold,
$\gamma$/hadron separation power and angular resolution of an ensemble of
detector units (Section 6). The tank candidates and dimensions presented here
are those that are most relevant and promising based on single tank and array
simulations.

##  2 Design Overview

A candidate Water Cherenkov Detector (WCD) design for SWGO with muon
separation potential is a double layered unit, as illustrated in Fig. 1, that
comprises the following building blocks:

Figure 1: Cylindrical double-layered WCD design comprising an upper chamber of
3.8 m diameter and 2.5 m depth, here with white walls and black bases (top and
bottom) and an entirely white lower chamber of 0.5 m depth. The upper chamber
contains an 8” PMT facing upwards, and the lower chamber has an 8” PMT facing
downwards. For illustration, a simulated muon (green track) is shown that
passes through both units and produces Cherenkov photons (red tracks).

  * •

Upper Chamber: A light-tight chamber with a lining that may be black or
reflective and a single upward-facing light sensor. This chamber provides
timing information and an estimate of total local particle energy per unit
area. The upward-facing light sensor ensures that non-reflected Cherenkov
photons with the smallest time dispersion are detected first.

  * •

Lower Chamber: A light-tight chamber with a highly reflective lining,
containing a single light sensor facing downwards (for improved uniformity of
response). This chamber enables muon tagging as only a small fraction of the
higher energy photons and electrons at ground level can punch through into the
lower chamber, while nearly all muons will pass through the entire detector
unit as the mean linear stopping power ($\rho\langle-dE/dx\rangle$) for muons
in water is only 2 MeV/cm.

The two PMTS are connected to each other with a PMT support such that one
faces upwards in the upper chamber, and the other faces downwards in the lower
chamber, as can be seen from Fig. 1. The reference design makes use of a 3.8 m
tank, motivated by the relative ease of road transportation of pre-fabricated
(e.g. rotomolded) units up to this size. The depths of the two chambers and
their reflective properties are investigated in Sections 4 and 5.

##  3 Overview of Simulation Tools

To simulate air showers, we use the CORSIKA 7.7400 simulation package Heck et
al. [11]. For the standard simulated event set, we select QGSJet-II.04
Ostapchenko [12] for energies above 80 GeV. UrQMD 1.3.1 Bass et al. [13],
Bleicher et al. [14] treats the low energy hadronic interactions and for
electromagnetic processes, we use the EGS4 electromagnetic model Nelson and
Namito [15].

We use GEANT4 Agostinelli et al. [16] within a simulation framework adapted
from that of the HAWC collaboration, to simulate the WCD response to the
secondary air shower particles from 20 m above the detector unit. This
simulation framework, called HAWCSim, has been extensively used for studies
related to HAWC and has been verified by the HAWC Collaboration Ahron S.
Barbera and Springer [17], Abeysekara and et al [18], Rehbein and et al [19].
The UNIFIED Levin and Moisan [20] model in GEANT4 is adapted to describe the
reflectivity of materials with a rough surface such as Tyvek (high
reflectivity - 92% at 450 nm, which from now on we refer to as ‘white’) Nozka
et al. [21] as used by the Pierre Auger Observatory Allekotte et al. [22] and
Polypropylene (low reflectivity - 10% at 450 nm, from now on referred to as
‘black’) as used by HAWC. The standard deviation of the distribution of the
micro-facet orientations, $\sigma_{\alpha}=0.17$ is set for the Tyvek surface.

A nominal water absorption length of 17 m at 400 nm is used in the
simulations. This water absorption length is reasonable when compared with the
numbers measured for the current detector arrays. A study from the LHAASO
collaboration using a custom developed apparatus reports attenuation lengths
around 16.8 m at 405 nm Li et al. [23]. Studies performed by the HAWC
collaboration shows attenuation lengths varying between 5 m and 16 m for 405
nm Abeysekara and et al [24].

We model the PMT in the simulations after the 8” HPK R5912-20 with a photo-
cathode quantum efficiency of 20% at 450 nm Hamamatsu [25]. The central 10”
PMT for one of the configurations in our comparisons (configuration B from
Section 4.1, similar to the HAWC central tanks) is modeled after the 10” R7081
with a photo-cathode efficiency of 30% Hamamatsu [25].

##  4 Upper Chamber Optimisation

For $\gamma$-ray induced extensive air showers (EAS) at typical detector
altitudes, the energy distributions of secondaries (in terms of number per log
energy interval) varies with particle type ($dN/dlogE$). This number density
peaks around $\sim$6 MeV for secondary photons, $\sim$20 MeV for electrons and
2–3 GeV for muons Schoorlemmer et al. [7]. A high detection probability for
particles of these energies is desirable for triggering and reconstructing
showers, combined with a precise determination of particle arrival time. Below
we discuss the impact of geometry and material properties on the performance
of the upper chamber, concerning particle detection efficiency (Section 4.1)
and arrival time measurement (Section 4.2). Chamber characteristics will also
influence the background trigger rates of each WCD unit (Section 4.3), that in
turn determine the trigger condition for the array and hence its energy
threshold.

###  4.1 Particle Detection Efficiency

The opening angle of the Cherenkov cone in water is $\sim$41°, and the pair
production length of high energy $\gamma$-rays is $\frac{9}{7}X_{0}$, where
$X_{0}$ is the radiation length, corresponding to $\sim$46 cm in water. The
depth of the chamber must be at least several radiation lengths to ensure that
all shower particles interact in the chamber. The diameter of the chamber then
determines the probability of collecting Cherenkov photons at the PMT, as
prompt photons or after some number of diffuse reflections. The WCD units of
both HAWC and LHAASO are lined with black material, with a water depth
${>}10\,X_{0}$ and a diameter-to-depth aspect ratio of 1.6–1.8.

Noting that prompt Cherenkov photons are important for timing (see Section
4.2) and that the chamber must be at least several radiation lengths
($\sim$4–5 $X_{0}\approx 1.5$–$1.8$ m) deep for calorimetric detection of
electromagnetic shower particles, we investigate the optimum chamber aspect
ratio. To illustrate the dependencies, we inject vertical 10 MeV $\gamma$-rays
uniformly across the top surface of a double-layered WCD unit with an 8” PMT
centered at the bottom. Fig. 2 shows the probability of detecting at least one
Cherenkov photon, as a function of chamber radius and depth, for chambers with
entirely white walls and with entirely black walls. As expected, the increased
photon path length and the nearly isotropic scattering of photons in the white
chamber increases the probability of light collection compared to a black
chamber. The white chamber provides detection efficiencies of 70% and more.
The detection probability decreases with increasing radius of the chamber. The
optimum depth varies with radius, between around 1 m to 2 m depth for the
white chamber of 1 m to 3 m radius and around 1.5 m to 2 m depth for the black
chamber.

The depth for optimum detection efficiency tends to be lower than the depth of
at least $5X_{0}$ required for good containment of electromagnetic shower
particles (red dashed line in Fig. 2), and the depth of at least 2.5 m
required for efficiently shielding the lower section from electromagnetic
particles (see Section 6.2). This shielding effect is illustrated in Fig. 3,
which shows the distribution of charge in the entirely white lower chamber,
for different upper chamber depths (3.8 m diameter and fixed 0.5 m lower
chamber depth) and different $\gamma$-ray energies. As the upper chamber depth
increases, the number of photo-electrons in the lower chamber is seen to
decrease.

Figure 2: Detection probability as a function of upper chamber radius and
depth, for vertical 10 MeV $\gamma$-rays injected across the top surface of
the chamber, for a chamber with entirely white walls (left) and a chamber with
entirely black walls (right). A dashed horizontal red line shows the depth
corresponding to $\sim$5 radiation lengths. Figure 3: The number of photo-
electrons in the lower chamber of a WCD unit with $\gamma$-ray energy ranging
from 10 MeV to 1 GeV for different upper chamber depths ($d_{u}$), again for
vertical 10 MeV $\gamma$-rays injected across the top surface of the chamber.
The mean number of photo-electrons ($\mu$) along with a dashed horizontal line
at 100 photo-electrons is shown for reference.The depth of the upper chamber
($d_{u}$) creates a shielding effect, as with increasing $d_{u}$ the number of
photo-electrons in the lower chamber decrease.

Furthermore, we investigate particle detection efficiency for four different
configurations of WCD units, detailed below.

  1. (A)

A white cylindrical double-layered WCD (3.8 m diameter and 2.5 m depth) with a
black top and an 8” PMT

  2. (B)

A cylindrical HAWC-like single-layered unit (7.3 m diameter and 4 m depth)
with black walls, a central 10” PMT and three 8” PMTs

  3. (C)

A LHAASO-like entirely black unit (5 m $\times$ 5 m square, 4.5 m depth) with
an open top and an 8” PMT111An 8” PMT is used for comparison and does not
reflect the actual PMT size(s) currently used in LHAASO WCD units

  4. (D)

A white cylindrical double-layered WCD unit (3.4 m diameter and 3.0 m depth)
with a black top and an 8” PMT

The chamber depth for configuration A is chosen based on the array
investigations from Section 6 (see Fig. 24) and is the reference configuration
for SWGO at the time of writing. Configuration D is a deeper and thinner
version of configuration A, as this variation is also being considered as a
viable alternative. The other two configurations are for options that would
essentially replicate the WCD units of the existing widefield observatories.
The central 10” PMT for configuration D is a higher quantum efficiency PMT
similar to HAWCAbeysekara et al. [1].

In order to investigate particle detection efficiency as a function of energy,
we inject vertical 5 MeV to 1 GeV $\gamma$-rays, uniformly distributed across
the top surface of these different configurations of WCDs. Fig. 4 shows a
comparison between the top chambers. The upper chamber of the double-layered
WCD (configurations A and D) has improved particle detection efficiency for
$\gamma$-rays over both HAWC-like and LHAASO WCD-like designs, due to the
reflective walls and – in case of the LHAASO-like design – the better ratio of
PMT area to chamber surface.

Figure 4: Detection probability for vertical 5 MeV to 1 GeV $\gamma$-rays
injected across the top surface of different WCD designs: (A) a cylindrical
white double-layered WCD (3.8 m diameter and 2.5 m depth) with a black top and
an 8” PMT; (B) a cylindrical single-layered HAWC-like unit (7.3 m diameter and
4 m depth) with black walls, a central 10” PMT and three 8” PMTs; (C) a LHAASO
WCD-like entirely black unit (5 m $\times$ 5 m square, 4.5 m depth) with an
open top and an 8” PMT; and (D) a white cylindrical double-layered WCD unit
(3.4 m diameter and 3.0 m depth) with a black top and an 8” PMT (note: An 8”
PMT is used here for comparison only and does not reflect the actual PMT(s)
size(s) currently used in LHAASO WCD units).

Multiple scattering of electrons pair produced by the incoming gammas, or in
later generations of a cascade, result in Cherenkov light emission in any
direction in the tank chamber. The reflectivity of the walls affects the
particle detection efficiency of the detector unit, as reflective walls result
in isotropisation of the emitted light. Entirely white walls will provide best
efficiency, but will result in Cherenkov photon arrival times with long tails,
the decay time scale governed by both the wall reflectivity and the water
transparency; at least the latter may vary over time and between detectors. We
consider combinations of white and black chamber walls in order to limit the
number of ‘late’ photons, such as a black top as used for the simulations
shown in Fig. 4. For the double-layered WCD (3.8 m diameter and 2.5 m depth,
option A from the different configurations), we compare all-white and all-
black chambers with chambers that have a black top, black bottom, or black
bases (i.e. top and bottom), for vertical 1 MeV to 1 GeV $\gamma$-rays
injected uniformly across the top surface, shown in Fig. 5. We see an
improvement in the particle detection efficiency of at least partially white
chambers over an entirely black chamber. Furthermore, as the number of photo-
electrons(pe) per unit energy scales like $\text{Area}_{\text{photo-
cathode}}/\text{Area}_{\text{chamber}}$ , we show in Tab. 1 the average number
of photo-electrons produced per $20$ MeV of $\gamma$-ray energy. Here a $20$
MeV $\gamma$-ray produces $\sim$1 pe in an entirely black chamber, while
producing ${\sim}3\times$ as many photo-electrons in an entirely white
chamber. The differences at such low energies for different material
combinations arise due to the effective production of a diffuse glow at the
top of the tank which a white top helps to deflect more towards the PMT. To
efficiently detect ${\leq}20$ MeV $\gamma$-rays would require reflective
materials and/or more photo-cathode efficiency/area.

Upper Chamber |  pe / 20 MeV
---|---
white |  2.9
white with black bases |  1.4
white with black bottom |  1.6
white with black top |  1.9
black |  1.0
Table 1: The average number of photo-electrons produced per 20 MeV
$\gamma$-ray energy in the upper chamber of a double-layered WCD unit (3.8 m
diameter and 2.5 m depth) with different materials. Figure 5: Detection
probability for vertical 1 MeV to 1 GeV $\gamma$-rays injected across the top
of the upper chamber of the double-layered WCD unit (A), for different chamber
wall materials.

In order to illustrate the performance of the calorimetric measurement of
shower particles, we show in Fig. 6 the relation between electromagnetic
energy and the number of photo-electrons and in Fig. 7 the resolution
$\sigma_{E}/E$ as a function of electromagnetic energy, for different
combinations of wall materials. As before, the all-white chamber behaves best
in terms of pe yield and energy resolution, the all-black chamber is worst,
while the mixed-wall chambers are intermediate. For incident particle energies
above $\sim$500 MeV, the lower chamber also becomes sensitive to
electromagnetic cascades, as shown in the same figure. Although the lower
chamber provides rather poor resolution (Fig. 7), it may help to extend the
dynamic range of the system, which is important for detector units close to
the shower core and/or in very high-energy showers.

Figure 6: The average number of photo-electrons as a function of
electromagnetic energy for upper chambers of double-layered WCD units (3.8 m
diameter and 2.5 m depth) with different materials. Also shown is the
corresponding value for the lower chamber (3.8 m diameter and 0.5 m depth) as
a dashed black line. Figure 7: Energy resolution as a function of
electromagnetic energy for upper chambers of double-layered WCD units (3.8 m
diameter and 2.5 m depth) with different materials. The dashed black line
shows the corresponding value for the lower chamber (3.8 m diameter and 0.5 m
depth).

The impact of the reflectivity of the upper chamber material on the angular
resolution and $\gamma$/hadron separation is discussed in later sections (cf.
6.2, 6.3), and is found to be modest.

###  4.2 Particle Arrival Time

The arrival time distributions of Cherenkov photons at the PMT determine the
time resolution that can (potentially) be achieved with a WCD unit; these time
distributions and the background rates (Section 4.3) also place important
requirements on the readout electronics of a unit, including the length of the
signal trace, buffering requirements, and trigger design. The aspect ratio and
the material choice of the WCD units determine these time distributions and
hence influence time resolution, which in turn impacts the achievable angular
resolution for showers as discussed in Section 6.3.

Figure 8: Timing distribution of photons arriving at the PMT for the four
different WCD configurations from Section 4.1. The photons are initiated by
vertical $\gamma$-rays at 80 MeV. Only the the timing distributions from the
central 10” PMT in the HAWC-like design and the PMT in the upper chamber of
the double-layered designs are shown.

In Fig. 8 we show Cherenkov photon arrival time distributions for the four
different WCD designs from the previous section, resulting from 80 MeV
vertical $\gamma$-rays. Reflective materials result in longer tails in the
timing distribution of both double-layered designs compared to the HAWC and
LHAASO-like designs which use entirely black walls.

Additionally as an example, we show in Fig. 9 the arrival time of Cherenkov
photon distributions resulting from 120 MeV vertical $\gamma$-rays injected
across the top surface of double-layered WCD units (3.8 m diameter and 2.5 m
depth, similar to configuration A) with a combination of different material
properties. The timing distributions show that a white walled chamber with one
or both of the top and bottom surfaces black reduces the tail of the arrival
distribution of photons by $\sim$30–40% at this energy.

Mostly prompt photons arrive at the PMT in the entirely black chamber while we
observe reflected light in chambers that have some combination of white
reflective material. Reflected light dominates the first arriving photons in
chambers with highly reflective material below $\sim$100 MeV. However, with
increased $\gamma$-ray energy, the first photons become dominated by direct
photons for all surface combinations (see Fig. 10).

Figure 9: Timing distribution of photons arriving at the upper chamber of a
double-layered WCD unit (3.8 m diameter and 2.5 m depth, similar to
configuration A) with different materials, for vertical $\gamma$-rays at 120
MeV. Shown also are the corresponding mean photon lifetimes.  Figure 10: The
fraction of events in which the first photon arriving at the PMT is a wall-
reflected photon, for the upper chamber of a double-layered WCD unit (3.8 m
diameter and 2.5 m depth) with different material properties and for vertical
$\gamma$-rays. The fraction is expressed relative to the all-black WCD.

Fig. 11 illustrates the impact of changing water
quality/absorption length on the arrival time distribution, for the case of a
3.8 m diameter, 2.5 m deep chamber with white walls and black top and bottom.
There is a very modest impact on the tail of the time distribution provided an
absorption length of $>$10 m can be maintained. A nominal absorption length of
17 m at 400 nm is used in the simulations.

Figure 11: Time distribution of photons arriving at the upper chamber of a
double-layered WCD unit (3.8 m diameter, 2.5 m depth, and white walls with
black bases) with water of varying absorption length, for vertical 120 MeV
$\gamma$-rays injected across the top surface. The mean photon lifetime in
each case is indicated in the figure legend.

###  4.3 Background WCD Trigger Rates

Air shower arrays typically trigger based on a coincidence of triggers from
individual units. The required minimal number of coincident units determines
the energy threshold of the array, and depends both on the trigger rate of
individual units and on the coincidence time window. The background trigger
rate of WCD units will depend on their detection probability for different
types of background particles, and hence on their geometry and wall materials;
we expect chambers with some white surfaces to experience increased background
rates with respect to entirely black chambers.

To investigate these rates, we use EXPACS/PARMA Sato [26] to calculate
terrestrial cosmic-ray fluxes and angular distributions of neutrons, protons,
electrons,
positrons, muons, anti-muons, and photons with energies ranging from
$1$–$10^{6}$ MeV at an altitude of 4900 m a.s.l. and location 13° 51’ 24” S,
71° 1’ 30” W (a location that well characterizes the magnetic effect at all
the sites under consideration for SWGO, see Doro et al. [27] for an overview).
We inject particles aimed at a hemisphere of 12 m radius centred at the
detector under test so that it can hit the detector from all sides (i.e 0 to
90 degree in zenith angle with a $sin\cdot cos$ distribution) to calculate the
trigger rates. Such a radius includes the edges of the surrounding detectors
to account for scattered particles and shielding effects. In order to include
shielding, a mini-array made up of 20 tanks with a separation of 0.12 m
between each tank was used. The expected single-pe trigger rates in the upper
chamber of double-layered WCD units (3.8 m diameter and 2.5 m depth) composed
of different materials are shown in Fig. 12, broken down into contributions
from different background particle species. Reflective walls increase the
detection probability of low energy particles resulting in an increased
background rate. Nevertheless, we expect an array of white-walled detectors to
provide the lowest energy threshold, once the coincidence level is adjusted to
obtain a negligible rate of array-level noise triggers. Note that these
background trigger rates do not include the single-pe thermal noise rate from
the PMT itself, or possible contributions from radioactive decays in the
water. These EXPACS/PARMA-based rates have been verified with CORSIKA air
shower simulations to generate ground-level particle rates. The rates at 4900
m are approximately ${\times}$1.5 the rates at an altitude of 4100 m, for the
same coordinates of an entirely white upper chamber of a double-layered WCD
unit (3.8 m diameter and 2.5 m depth), see Fig. 13.

Figure 12: Noise trigger rates at single pe threshold for the upper chamber of
double-layered WCD unit (3.8 m diameter and 2.5 m depth) comprising an 8” PMT
and various combination of materials (b = black, w = white), at an altitude of
4900 m. Contributions from $e^{\pm}$, $\mu^{\pm}$, $\gamma$, $p^{+}$ and $n$
are shown.

The same method was used to predict the noise trigger rates of HAWC tanks,
specifically of the central 10” PMT in a tank in the HAWC central array and an
HAWC outrigger tank, at the actual HAWC altitude of 4100 m. The rates, shown
in Figure 14 along with an afterpulsing rate prediction, agree with the
measured HAWC rates, even though contributions from dark rates that varies
around 1-3 kHz with temperature, voltage and PMT dependenceHamamatsu [25] are
not included in the method. The afterpulsing rates were estimated with a
first-order calculation that assumes an afterpulse probability of 0.02% for
the 10” PMTs, as reported for these PMTs for a 2013 study for the Double Chooz
experiment Haser et al. [28]. The prediction can vary, as PMTs can have higher
afterpulsing probabilities of 10% if they are degraded, which would mean that
the afterpulse rate is higher than the shown amount for some of the tanks.
Indeed, the HAWC data shows a large spread. The central 10” PMTs of HAWC were
reported to have a hit rate of 40-50 kHz Abeysekara et al. [1]. The outriggers
of HAWC were reported to have rates around 4-8 kHz Marandon et al. [29], which
agrees with our predictions with this method.

Figure 13: Noise trigger rates at single pe threshold for an entirely white
upper chamber of double-layered WCD unit (3.8 m diameter and 2.5 m depth)
comprising an 8” PMT at various altitudes. Contributions from $e^{\pm}$,
$\mu^{\pm}$, $\gamma$, $p^{+}$ and $n$ are shown. Figure 14:  Noise trigger
rates at single pe threshold for two types of HAWC tanks obtained with HAWCSim
at 4100 m. Contributions from $e^{\pm}$, $\mu^{\pm}$, $\gamma$, $p^{+}$ and
$n$ are shown _Left:_ Rates for the central PMT of a tank in the central HAWC
array. _Right:_ Rates for a HAWC outrigger tank.

##  5 Lower Chamber Optimisation

In the dual-layer approach muons are identified based on the signal from the
lower chamber, or, more generally, by comparing the signals in the upper and
lower chambers. Since the lower chamber is not used for timing, we assume
white-walled lower chambers for optimal light yield. The depth of the lower
chamber, combined with the photosensor area, determines the light yield. While
in principle a few detected photons per muon are sufficient to tag muons with
reasonable efficiency, our studies of algorithms for muon identification in
showers (Section 6.2) indicate that larger signals are desirable to reduce the
number of misidentified muons. Here, we discuss how the muon signal depends on
the geometry of the lower chamber.

Figure 15: The number of photo-electrons due to incident $\mu^{-}$ of 100 MeV
to 2 GeV in a WCD unit with a white-walled (3.8 m diameter) lower cell of
various depths. Dashed lines show 32 pe and 1 GeV.

###  5.1 Depth and Reflectivity of the Lower Chamber

We study lower chambers with white walls and with a PMT placed at the centre
of the top of the chamber facing downwards. The PMT base protrudes into the
upper chamber so that only the active cathode area of the PMT is visible in
the lower chamber, for best photon collection efficiency and hence minimal
chamber depth. To investigate the optimum depth for the lower chamber,
vertical muons with energies from 100 MeV to 2 GeV are injected uniformly
across the surface of the 3.8 m diameter upper chamber. The lower chamber has
the same diameter to allow partitioning of a single cylindrical detector unit.

Fig. 15 shows the photo-electrons collected as a function of energy for
varying depths of the lower chamber with an 8” PMT. There is a uniform
response for muons above $\sim$600 MeV as expected since muons lose $\sim$2
MeV/cm in water. Figure 16 shows that below a depth of $\sim$0.5 m, the muon
signal has a much larger spread compared to depths including and above
$\sim$0.5 m. Moreover, above a depth of $\sim$0.5 m, the muon signal is
reliably $>$10 photo-electrons for $>$1 GeV muons, and the signal increases
roughly proportionally to the track length of the impinging particle as
expected (see Fig. 16). It is seen that a depth of $\sim$0.5 m and above would
give reliable muon signals, and since a lower chamber that is smaller in depth
would minimize costs (although most of the cost is driven by the upper
chamber), the studies here suggest that the depth for the lower chamber should
be at least $\sim$0.5 m.

Figure 16: The number of photo-electrons in a WCD unit with an entirely white
lower chamber, for 1 GeV $\mu^{-}$ energy and various chamber depths. Figure
17: The number of photo-electrons in a WCD unit with an entirely white upper
(3.8 m diameter and 2..5 m depth) and lower (3.8 m diameter and 0.5 m depth)
chamber. On the left, for $\mu^{-}$ with energy ranging from 100 MeV to 2 GeV;
on the right, for $\gamma$-rays with energy ranging from 1 MeV to 2 GeV.
Dashed lines show 80 pe and 1 GeV in the upper plots and 32 pe and 1 GeV in
the lower plots.

To explore the separation power of such a unit, we compare vertical
$\gamma$-rays with energies from 1 MeV to 2 GeV in an entirely white WCD unit
(3.8 m diameter, 2.5 m upper and 0.5 m lower depth) with vertical $\mu^{-}$
with energies from 100 MeV to 2 GeV in an identical WCD unit (see Fig. 17).
Muon identification is possible as the number of photo-electrons detected in
the two chambers is constant above some energy threshold and remains fairly
stable until very high energies where effects such as bremsstrahlung will need
to be considered. Additionally, the peak of the muon energy distribution in
air showers is around 2–3 GeV Schoorlemmer et al. [7]. The ratio of photo-
electrons in the two chambers will enable muon selection on a tank-by-tank
basis. Other shower information such as the location of the shower core
relative to the detector unit and the number of units hit is also beneficial
for muon identification.

##  6 Array Simulations

Figure 18: Sketch of the simulated reference configuration for SWGO comprising
cylindrical double-layered WCDs with a dense inner array ($\sim$80%) and
sparser outer array ($\sim$8%) The inner array has 0.12 m separation between
the individual units.

To relate the performance of individual WCD units to the performance of an
array as a whole and in particular to the achievable angular resolution and
$\gamma$/hadron separation power, we carried out simulations based on the SWGO
reference design Schoorlemmer et al. [30]. We continue using the _HAWCSim_
tool that has been validated by the HAWC collaboration, as explained in
Section 3. We simulate an inner array with a high fill factor of $\sim$80%
spread over $\approx$80,000 m2 and a sparser outer array with a fill factor of
8% spread over $\approx$221,000 m2, as shown in Fig. 18. The array is placed
at high elevation (4900 m.a.s.l). The dense inner array serves to enhance the
sensitivity for low to mid energy $\gamma$-rays and also to increase the muon
sensitive area. The sparser outer array is designed for a large collection
area at the highest energies, but will not be used in the analysis presented
here.

###  6.1 Estimating Effective Area

To be able to compare the $\gamma$-ray induced air shower effective area for
different material combinations, we first simulate proton (background) induced
showers with a spectral index of -2, an energy range from 0.001–30 TeV and
zenith angle from $0^{\circ}\;{\rm to}\;60^{\circ}$, at an altitude of 4900
m.a.s.l. and shower core scattered over a radius of 2.5 km. These simulations
are then weighted by the cosmic-ray flux Lipari and Vernetto [31] to reproduce
the correct spectrum. As the core range is large compared to the detector
size, we can further integrate to obtain approximate array trigger rates at
different hit thresholds (see Fig. 19).

Figure 19: Background array rates for an array of double-layered WCDs (3.8 m
diameter, 2.5 m upper depth, and 8” PMT) composed of different materials.

We then simulate $\gamma$-ray induced air showers with a spectral index of -2
for energy range from 0.1–5 TeV, zenith angle from $0^{\circ}\;{\rm
to}\;60^{\circ}$, at an altitude of 4900 m.a.s.l. and shower core scattered
over a radius of 2.5 km. From the proton simulations, the required threshold
in the number of array hits is derived, such that the array trigger rate is
100 kHz (ensuring at least a few 10s of tanks trigger). With this threshold,
the $\gamma$-ray simulations are used to derive the gamma-ray effective area
as a function of energy (see Fig. 20). The simulations show that double-
layered WCDs (3.8 m diameter, 2.5 m upper depth, and 8” PMT) with some
reflective surfaces for the upper chamber provide slightly higher effective
areas at low energies, compared to entirely black upper chambers. It should be
noted however that more statistics, including hit timing information and large
core range would be necessary to obtain realistic rates.

Figure 20: $\gamma$-ray induced air shower effective area for an array of
double-layered WCDs (3.8 m diameter, 2.5 m upper depth, and 8” PMT) composed
of different materials for the upper chamber.

###  6.2 Identifying Hadronic Showers

The double-layered WCD design provides an alternate way of vetoing hadron-
induced showers compared to single-chambered WCD designs. The critical
parameter that we will explore in the following is the depth of the upper
chamber, that determines the level of electromagnetic punch-through into the
lower chamber.

For this purpose a template-based maximum log-likelihood method was
implemented to discriminate between $\gamma$-ray and hadron-induced air
showers. Given a known core location and air shower direction, templates for
the charge in the upper and lower chambers in an array of double-layered WCDs,
similar to the distribution shown in Fig. 21 are generated for vertical air
showers. Separate templates are generated for $\mu^{\pm}$ and for other
charged particles. Next, we test simulation events, where we assign a
likelihood value based on the charge deposited in the two chambers for each
secondary particle impinging on an individual tank. Subsequently, we calculate
a Likelihood Ratio (LR) to tag those particles with a $LR>0$ as candidate
muons and compare them to the Monte-Carlo truth as shown in Fig. 22. The
figure shows that muon identification is effective for muon events with high
upper chamber signal.

Figure 21: Shower templates for electromagnetic and muon energy deposition,
showing charge in the upper and lower chambers of entirely white double-
layered WCDs for proton showers, and excluding the region within 20 m of the
shower core. Figure 22: Performance of the likelihood ratio based
identification of muons, illustrated as a function of signal amplitude in the
upper chamber for vertical $\gamma$ and proton-initiated showers on an array
of entirely white double-layered WCDs. Identification is very effective for
those muon events in which the upper chamber signal is sufficiently high.

This muon tagging ability directly translates into $\gamma$/hadron separation
efficiency, given the relative abundance of muons in hadron-initiated showers.
Using the number of detector units hit as a proxy for the shower’s energy, we
can subsequently distinguish $\gamma$-ray and hadron-initiated showers from
the difference in the number of identified muons, as shown in Fig. 23 for
vertical 2 TeV $\gamma$-ray and 5 TeV proton induced showers. The proton
shower energy is chosen such the the average number of WCD hits is very
similar for the two sets of events. To reduce misidentification due to punch-
through of electromagnetic particles from the upper to lower chamber close to
the shower core, we exclude WCDs within 20 m from the core.

Figure 23: Probability of detecting muons in detector units more than 20 m
from the shower core located at the centre of the array, for $\gamma$-ray- and
hadron-initiated showers detected with entirely white double-layered WCDs in
the dense inner array of the SWGO Reference Configuration. Here a black dashed
line is shown where the SNR
($\epsilon_{\gamma}/\sqrt{\epsilon_{\text{proton}}}$) is maximum for
reference.

From the known core distance and air shower direction, we generate templates
of charge in the upper and lower chambers in an array of entirely white
double-layered WCDs (3.8 m diameter and 2.5 m upper chamber depth and 0.7 m
lower chamber depth) for vertical 5 TeV proton-induced air showers.

After tagging different particle species for 2 TeV $\gamma$-ray and 5 TeV
proton induced showers, we find a $\gamma$/hadron separation efficiency
varying the upper chamber depth. The signal-to-noise-ratio (SNR) is
$\epsilon_{\gamma}/\sqrt{\epsilon_{\text{proton}}}$, where $\epsilon$ is
efficiency. Above an upper chamber depth of 2.5 m, the SNR plateaus for air
showers at 0∘ zenith angle. Inclined showers are also expected to plateau, but
as individual secondary particles from inclined showers can penetrate through
the sides and multiple tanks, there is a decrease in overall SNR. Including
the neighbouring tanks would mitigate the decrease in SNR for such inclined
showers. Increasing detector unit radii or increasing the fill factor could
also be beneficial for shielding from side-penetrating particles; one could
also consider filling up the space between tanks with ground material as
absorber. However, a complete study of inclined showers is beyond the scope of
this paper.

Attempting to account for the cost scalings for tanks with depth, volume of
water, and the fixed cost of photosensors, we find that for vertical showers
the SNR-to-cost ratio peaks at $\sim$2.5 m, and this depth is relatively
insensitive to the cost assumptions made (see Fig. 24). Based on these
findings, the upper chamber depth of 2.5 m has been used in the comparisons of
Section 4 and while investigating different chamber materials.

Figure 24: Ratio of $\gamma$/hadron separation SNR to nominal cost in an array
with entirely white double-layered WCDs with varying upper chamber depth,
given different assumptions of the relative cost of water to photosensors and
tanks, assuming a fixed number of tanks. Each WCD unit has a diameter of 3.8 m
and a fixed lower chamber depth of 0.7 m.

To show the effect of various materials in the upper chamber for double-
layered WCDs (3.8 m diameter, 2.5 m upper depth, 0.5 m lower depth, with 8”
PMTs), we again implement a template-based maximum likelihood method to
discriminate between $\gamma$-ray and hadron induced air showers for the array
shown in Fig. 18. We simulate an ensemble of vertical $\gamma$-ray and proton
initiated showers following an E-2 spectrum up to 100 TeV energy, with shower
core location at the centre of the array. Fig. 25 shows the derived
$\gamma$/hadron separation efficiency (excluding a 40 m region around the
core) for different material combinations. It is seen that a background
rejection power of $\sim$103 can be achieved at reconstructed energies (a
simple model taking into account only the total number of photo-electrons seen
in all of the upper chambers of the array) of a few TeV (2- 7 TeV) with high
$\gamma$-ray efficiency. The separation power improves in all cases with one
or more white surface(s).

Figure 25: $\gamma$/hadron separation efficiency for an array of double-
layered WCDs (3.8 m diameter and 2.5 m depth) with $\sim$80% fill factor
varying material reflectivity, using an exclusion region of 40 m around the
core location for reconstructed Gamma showers of 2 to 7 TeV.

###  6.3 Angular Reconstruction

To reconstruct the direction of the air shower, we use the time of the first
photon in each upper WCD without applying transition time spread in the PMT,
or electronics time resolution. We first obtain a guess of the direction from
a least-squares algorithm using the time difference of arrival of the photons
between detector units.

The final best-fit direction is obtained via a likelihood fit to the arrival
times in each hit detector unit. Parameterized distribution functions describe
the distribution in energy, core distance and arrival time of shower particles
at ground. PDFs for the time distribution are obtained by fitting Landau
distributions to the arrival times binned as function of distance to the
shower core and total charge, to obtain mean and width parameters for each bin
Hofmann [32]. A shower arrival time, $t_{0}$, is defined as the absolute
arrival time of the first electromagnetic particle in the plane perpendicular
to the shower axis. To separate the effects of the direction (i.e. timing) fit
and the core position fit, and to minimise additional complexity due to array
edges, we generate vertical $\gamma$-rays impacting at the centre of the
array, and assume the core position to be known (as in practice the core
location can be very precisely determined for well-contained events). A three-
parameter (time offset and direction) likelihood fit is performed (using
MINUIT James [33]) to obtain the reconstructed shower directions and hence
determine the angular resolution (see Fig. 26). As expected, the angular
resolution improves with increasing energy. Differences between different wall
options are very modest. The additional late-arriving Cherenkov photons that
are a feature of white-walled chambers do not deteriorate the angular
resolution as they are properly accounted for in the likelihood functions used
in the fit; these likelihood functions are adapted for each configuration. In
all cases, the plot suggests that the angular resolution obtained is superior
to the LHAASO and HAWC observatory angular resolutions Aharonian et al. [2],
Abeysekara et al. [1]. Whilst this result represents a very idealised case,
with showers from zenith landing at the centre of the array, with perfect
timing resolution and no noise, it nonetheless illustrates that such arrays
can potentially achieve an angular resolution much better than existing
instruments of this type, regardless of the surface reflectivity of individual
detector units. Comparisons with different fill factors and reconstruction
algorithms are beyond the scope of this paper.

Figure 26: Angular resolution for 1–100 TeV vertical $\gamma$-rays at the
centre of an array of double-layered WCDs (3.8 m diameter and 2.5 m depth)
with $\sim$80% fill factor – for a range of options on surface reflectivity.

##  7 Summary & Outlook

We have studied a double-layered WCD array, aiming to improve both the energy
threshold and the $\gamma$-ray/hadron separation compared to LHAASO- and
HAWC-like designs. Each detector unit comprises two chambers with black or
diffuse reflective wall linings and a PMT in each chamber. The upper PMT
facing upwards is intended for timing and energy determination, while the
lower PMT facing downwards will enable muon tagging.

Investigating different options for chamber aspect ratio (depth-to-diameter)
and reflectivity shows that a compact and highly reflective upper chamber
lowers the energy threshold for the impinging particles and, in conjunction
with a lower chamber, benefits $\gamma$/hadron separation. The double-layer
design is promising in achieving background rejection power of $\sim$103 at
energies of a few TeV. There is an increase in the tail in the timing
distributions with one or more white surface in the upper chamber, however,
these late-arriving Cherenkov photons do not deteriorate the angular
resolution obtained. In the idealized case of vertical showers at the centre
of the array an angular resolution of several arc-minutes at 10 TeV was
obtained. A partially reflective upper chamber is strongly motivated by these
benefits in threshold and background rejection, with no negative impact on
resolution. As for the chamber depths, an upper chamber depth of $\sim$2.5 m
maximizes the SNR-to-cost ratio and provides sufficient shielding for the
lower section, meanwhile a lower chamber depth of $\sim$0.5 m is the minimum
required for reliable muon signals.

Although the final optimization of a unit detector for SWGO requires finalized
analysis algorithms and life cycle cost estimates that are beyond the scope of
this paper, these studies show that a $\sim$3 m deep, $\sim$4 m diameter
double-layered tank with some reflective material in the upper chamber is a
promising option in terms of performance and cost-effectiveness.

##  8 Acknowledgements

We thank our colleagues within the SWGO collaboration for many helpful
discussions and the use of the common shared software framework. We thank the
HAWC collaboration for providing the AERIE software used here to SWGO.

## References

  * Abeysekara et al. [2017] A. U. Abeysekara, et al.,  Observation of the crab nebula with the hawc gamma-ray observatory,  The Astrophysical Journal 843 (2017) 39. URL: https://doi.org/10.3847/1538-4357/aa7555. doi:10.3847/1538-4357/aa7555.
  * Aharonian et al. [2021] F. Aharonian, et al.,  Observation of the crab nebula with lhaaso km2a a performance study,  Chinese Physics C 45 (2021) 025002. URL: http://dx.doi.org/10.1088/1674-1137/abd01b. doi:10.1088/1674-1137/abd01b.
  * Abreu et al. [2019] P. Abreu, et al., The southern wide-field gamma-ray observatory (swgo): A next-generation ground-based survey instrument for vhe gamma-ray astronomy, 2019. arXiv:1907.07737.
  * Hinton [2021] J. Hinton,  The southern wide-field gamma-ray observatory: Status and prospects,  Proceedings of 37th International Cosmic Ray Conference PoS(ICRC2021) (2021).
  * Conceição, R. et al. [2021] Conceição, R., González, B. S., Guillén, A., Pimenta, M., Tomé, B.,  Muon identification in a compact single-layered water cherenkov detector and gamma/hadron discrimination using machine learning techniques,  Eur. Phys. J. C 81 (2021) 542. URL: https://doi.org/10.1140/epjc/s10052-021-09312-4. doi:10.1140/epjc/s10052-021-09312-4.
  * Sinnis [2007] G. Sinnis,  High altitude water water cherenkov telescope,  in: 1st Roma International Conf. on Astro-Particle Physics (RICAP 07), Rome, Italy, 2007\.
  * Schoorlemmer et al. [2019] H. Schoorlemmer, J. Hinton, R. López-Coto,  Characteristics of extensive air showers around the energy threshold for ground-particle-based $$$\backslash$gamma $$-ray observatories,  The European Physical Journal C 79 (2019) 427. URL: https://doi.org/10.1140/epjc/s10052-019-6942-x. doi:10.1140/epjc/s10052-019-6942-x.
  * Assis et al. [2017] P. Assis, et al.,  LATTES: a new gamma-ray detector concept for South America,  EPJ Web Conf. 136 (2017) 03013. doi:10.1051/epjconf/e2016-62145-8. arXiv:1703.09254.
  * Atkins and et al [2003] R. Atkins, et al,  Observation of TeV Gamma Rays from the Crab Nebula with Milagro Using a New Background Rejection Technique,  apj 595 (2003) 803–811. doi:10.1086/377498. arXiv:astro-ph/0305308.
  * Letessier-Selvon et al. [2014] A. Letessier-Selvon, P. Billoir, M. Blanco, I. C. Mariş, M. Settimo,  Layered water cherenkov detector for the study of ultra high energy cosmic rays,  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 767 (2014) 41–49. URL: https://www.sciencedirect.com/science/article/pii/S0168900214009589. doi:https://doi.org/10.1016/j.nima.2014.08.029.
  * Heck et al. [1998] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz, T. Thouw, CORSIKA: A Monte Carlo code to simulate extensive air showers, 1998.
  * Ostapchenko [2011] S. Ostapchenko,  Monte carlo treatment of hadronic interactions in enhanced pomeron scheme: Qgsjet-ii model,  Phys. Rev. D 83 (2011) 014018. URL: https://link.aps.org/doi/10.1103/PhysRevD.83.014018. doi:10.1103/PhysRevD.83.014018.
  * Bass et al. [1998] S. Bass, M. Belkacem, M. Bleicher, M. Brandstetter, L. Bravina, C. Ernst, L. Gerland, M. Hofmann, S. Hofmann, J. Konopka, G. Mao, L. Neise, S. Soff, C. Spieles, H. Weber, L. Winckelmann, H. Stöcker, W. Greiner, C. Hartnack, J. Aichelin, N. Amelin,  Microscopic models for ultrarelativistic heavy ion collisions,  Progress in Particle and Nuclear Physics 41 (1998) 255–369. URL: https://www.sciencedirect.com/science/article/pii/S0146641098000581. doi:https://doi.org/10.1016/S0146-6410(98)00058-1.
  * Bleicher et al. [1999] M. Bleicher, E. Zabrodin, C. Spieles, S. A. Bass, C. Ernst, S. Soff, L. Bravina, M. Belkacem, H. Weber, H. Stöcker, et al.,  Relativistic hadron-hadron collisions in the ultra-relativistic quantum molecular dynamics model,  Journal of Physics G: Nuclear and Particle Physics 25 (1999) 1859–1896. URL: http://dx.doi.org/10.1088/0954-3899/25/9/308. doi:10.1088/0954-3899/25/9/308.
  * Nelson and Namito [1990] W. R. Nelson, Y. Namito, The EGS4 Code System: Solution of gamma-ray and electron transport problems, Nuclear Energy Data Center, Japan, 1990.
  * Agostinelli et al. [2003] S. Agostinelli, et al. (GEANT4),  GEANT4–a simulation toolkit,  Nucl. Instrum. Meth. A 506 (2003) 250–303. doi:10.1016/S0168-9002(03)01368-8.
  * Ahron S. Barbera and Springer [2017] D. B. K. Ahron S. Barbera, R. W. Springer,  Simulation of near horizontal muons and muon bundles for the hawc observatory with corsika,  2017\.
  * Abeysekara and et al [2015] A. Abeysekara, et al,  Vamos: A pathfinder for the hawc gamma-ray observatory,  Astroparticle Physics 62 (2015) 125–133. URL: https://www.sciencedirect.com/science/article/pii/S0927650514001285. doi:https://doi.org/10.1016/j.astropartphys.2014.08.004.
  * Rehbein and et al [2017] F. Rehbein, et al,  Preliminary cosmic ray results from the hawc’s eye telescopes,  2017\.
  * Levin and Moisan [1996] A. Levin, C. Moisan,  A more physical approach to model the surface treatment of scintillation counters and its implementation into detect,  in: 1996 IEEE Nuclear Science Symposium. Conference Record, volume 2, 1996, pp. 702–706 vol.2. doi:10.1109/NSSMIC.1996.591410.
  * Nozka et al. [2011] L. Nozka, M. Pech, H. Hiklova, D. Mandat, M. Hrabovsky, P. Schovanek, M. Palatka,  Brdf profile of tyvek and its implementation in the geant4 simulation toolkit,  Opt. Express 19 (2011) 4199–4209. doi:10.1364/OE.19.004199.
  * Allekotte et al. [2008] I. Allekotte, A. Barbosa, P. Bauleo, C. Bonifazi, B. Civit, C. Escobar, B. García, G. Guedes, M. Gómez Berisso, J. Harton, et al.,  The surface detector system of the pierre auger observatory,  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 586 (2008) 409–420. URL: http://dx.doi.org/10.1016/j.nima.2007.12.016. doi:10.1016/j.nima.2007.12.016.
  * Li et al. [2018] C. Li, G. Xiao, S. Feng, L. Wang, X. Li, X. Zuo, N. Cheng, H. Wang, B. Gao, Z. Duan, J. Liu, H. He, M. Saeed, LHAASO Collaboration,  An apparatus to measure water optical attenuation length for LHAASO-MD,  Nuclear Instruments and Methods in Physics Research A 892 (2018) 122–126. doi:10.1016/j.nima.2018.03.029.
  * Abeysekara and et al [2015] A. U. Abeysekara, et al,  SEARCH FOR GAMMA-RAYS FROM THE UNUSUALLY BRIGHT GRB 130427a WITH THE HAWC GAMMA-RAY OBSERVATORY,  The Astrophysical Journal 800 (2015) 78. URL: https://doi.org/10.1088/0004-637x/800/2/78. doi:10.1088/0004-637x/800/2/78.
  * Hamamatsu [2019] Hamamatsu, Large Photocathode Area Photomuliplier Tubes, 2019. URL: https://www.hamamatsu.com/resources/pdf/etd/LARGE_AREA_PMT_TPMH1376E.pdf.
  * Sato [2016] T. Sato,  Analytical model for estimating the zenith angle dependence of terrestrial cosmic ray fluxes,  PLOS ONE 11 (2016) 1–22. URL: https://doi.org/10.1371/journal.pone.0160390. doi:10.1371/journal.pone.0160390.
  * Doro et al. [2022] M. Doro, A. Moraes, M. Santander, M. Giammarchi, D. Mandat, L. M. Mendes, J. Bellido, C. Dib, S. May, A. Reisenegger, A. C. Rovero, A. Sandoval, F. Schüssler, R. W. Springer, I. Torres, J. Vicha, SWGO Collaboration,  The search for high altitude sites in South America for the SWGO detector,  in: 37th International Cosmic Ray Conference. 12-23 July 2021. Berlin, 2022, p. 689.
  * Haser et al. [2013] J. Haser, F. Kaether, C. Langbrandtner, M. Lindner, S. Lucht, S. Roth, M. Schumann, A. Stahl, A. Stüken, C. Wiebusch,  Afterpulse measurements of R7081 photomultipliers for the Double Chooz experiment,  Journal of Instrumentation 8 (2013) P04029. doi:10.1088/1748-0221/8/04/P04029. arXiv:1301.2508.
  * Marandon et al. [2019] V. Marandon, A. Jardin-Blicq, H. Schoorlemmer,  Latest news from the HAWC outrigger array,  volume 36 of International Cosmic Ray Conference, 2019, p. 736. arXiv:1908.07634.
  * Schoorlemmer et al. [2021] H. Schoorlemmer, et al.,  Simulating the performance of the southern wide-view gamma-ray observatory,  in: Proceedings of Science (ICRC 2021), 903, Berlin, Germany, 2021\.
  * Lipari and Vernetto [2019] P. Lipari, S. Vernetto,  The shape of the cosmic ray proton spectrum,  arxiv: High Energy Astrophysical Phenomena (2019).
  * Hofmann [2020] W. Hofmann,  On angular resolution limits for air shower arrays,  Astropart. Phys. 123 (2020). doi:10.1016/j.astropartphys.2020.102479.
  * James [1994] F. James, MINUIT Function Minimization and Error Analysis: Reference Manual Version 94.1, 1994.