# Enhanced Anomaly Detection in Automotive Systems Using SAAD: Statistical
Aggregated Anomaly Detection

1st Dacian Goina  Faculty of Mathematics
and Computer Science
West University of Timisoara
Timisoara, Timis
<EMAIL_ADDRESS>2nd Eduard Hogea  Faculty of Mathematics
and Computer Science
West University of Timisoara
Timisoara, Timis
<EMAIL_ADDRESS>3rd George Maties  Faculty of Electrical Engineering
Gheorghe Asachi
Technical University of Iasi
Iasi, Iasi
<EMAIL_ADDRESS>
###### Abstract

This paper presents a novel anomaly detection methodology termed Statistical
Aggregated Anomaly Detection (SAAD). The SAAD approach integrates advanced
statistical techniques with machine learning, and its efficacy is demonstrated
through validation on real sensor data from a Hardware-in-the-Loop (HIL)
environment within the automotive domain. The key innovation of SAAD lies in
its ability to significantly enhance the accuracy and robustness of anomaly
detection when combined with Fully Connected Networks (FCNs) augmented by
dropout layers. Comprehensive experimental evaluations indicate that the
standalone statistical method achieves an accuracy of 72.1%, whereas the deep
learning model alone attains an accuracy of 71.5%. In contrast, the aggregated
method achieves a superior accuracy of 88.3% and an F1 score of 0.921, thereby
outperforming the individual models. These results underscore the
effectiveness of SAAD, demonstrating its potential for broad application in
various domains, including automotive systems.

###### Index Terms:

anomaly detection, HIL, machine learning, automotive systems, fully connected
networks, statistical methods

##  I Introduction

Figure 1: Architecture of the SAAD(Statistical Aggregated Anomaly Detection)
framework, showcasing the integration of statistical anomaly detection with a
Fully Connected Network (FCN) for enhanced anomaly detection in automotive
systems. The architecture includes data preprocessing, artificial labeling,
deep learning model training, and result aggregation to achieve improved
accuracy and robustness. The effectiveness of SAAD is demonstrated with two
samples, X and Y. The FCN initially predicts that both are anomalies, with
sample X having high confidence and sample Y having low confidence. The
statistical method, using $k=2$, detects sample X as an anomaly and sample Y
as not. Aggregation is applied to obtain the final answer, depending on the
outputs of the models and the confidence of the deep learning one, combining
the strengths of both methods.

Anomaly detection is a critical task in various fields, particularly in the
automotive industry, where it plays a significant role in ensuring the
reliability and safety of systems. With the advent of advanced machine
learning (ML) techniques, anomaly detection has seen substantial improvements
in accuracy and efficiency. Machine learning models, especially deep learning
architectures, have shown promising results in identifying abnormal patterns
that could indicate potential faults or cyber-attacks in automotive systems
[1, 2].

In the context of automotive systems, anomaly detection is crucial for
maintaining the operational integrity of both traditional and autonomous
vehicles. Autonomous vehicles, in particular, rely heavily on a multitude of
sensors to navigate and operate safely. These sensors provide continuous
streams of data that need to be monitored in real-time for any signs of
anomalies that could compromise the vehicle’s safety and functionality [3, 4].
Machine learning models have been employed to analyze these data streams,
allowing for the early detection of issues and the implementation of
corrective actions before any adverse events occur [5].

To validate and verify the performance of control systems under real-time
conditions, Hardware-in-the-Loop (HIL) testing is an essential technique used
in the industry. HIL testing involves integrating actual hardware components
with simulated environments, providing a highly controlled and repeatable
setting to assess the behavior of automotive systems. This method is
particularly beneficial for anomaly detection, as it allows for the collection
of high-fidelity data that can be used to train and test machine learning
models under various scenarios [6]. Additionally, HIL testing helps mitigate
the common issue of imbalanced data in anomaly detection datasets by enabling
the generation of synthetic data to balance the dataset.

Evaluating anomaly detection algorithms often requires specific metrics
tailored for time series data, as classical precision and recall metrics may
not capture the nuanced performance characteristics required for detecting
anomalies in continuous data streams [7]. Advanced evaluation metrics consider
the duration and timing of anomalies, providing a more robust assessment of an
algorithm’s performance in real-world conditions.

Recent studies have explored various machine learning approaches for anomaly
detection in different domains. For instance, [8] utilized explainable random
forests to detect anomalies in oil wells, achieving high accuracy and
demonstrating the effectiveness of explainable AI in critical applications.
Another study [9] proposed a clustered deep one-class classification model for
detecting anomalies in cyber-physical systems, showing significant performance
improvements over traditional models. A survey on anomaly detection in IoT
time-series data highlighted the unique challenges and proposed solutions for
implementing these techniques in rapidly expanding fields like the Internet of
Things [10]. Additionally, a comprehensive review on outlier detection in time
series data provided a structured overview of unsupervised techniques and
their applications [11].

The integration of ML in the industry has also yielded impressive results in
structural health monitoring [12, 13]. Building on these advancements, recent
studies have expanded the field towards hybrid-assisted models that leverage
logical rules. These neuro-symbolic approaches combine deep learning with
knowledge-based systems to enhance performance. Research has shown that these
hybrid models, which incorporate domain-specific knowledge, significantly
outperform standard deep learning models by providing more accurate and
explainable results [14, 15]. Additionally, innovative methods such as
GraphLIME have demonstrated significant improvements in model interpretability
and performance in predictive tasks, further showcasing the potential of
combining explainable AI techniques with deep learning [16]. The
interpretability aspect is highly valuable, as seen in [17], where the
insights were used to understand and solve the feature distribution issue that
occurred with the generated pseudo-features. A comprehensive comparison in
[18] highlights both the advantages of hybrid models and the challenges of
creating and scaling them. These hybrid approaches, which synergize
complementary components, form the cornerstone of the conceptual framework
presented in this paper.

In this paper, we introduce Statistical Aggregated Anomaly Detection (SAAD), a
novel anomaly detection method designed to effectively integrate statistical
techniques for outlier detection with the performance of supervised machine
learning models. Our approach specifically aims to enhance the accuracy and
robustness of Fully Connected Networks (FCNs) in anomaly detection tasks.
Utilizing real-world data collected from Hardware-in-the-Loop (HIL)
environments, we demonstrate the effectiveness of SAAD in identifying
anomalies with high precision and recall. Additionally, we illustrate how our
statistical method can augment the performance of FCNs, rendering them more
reliable for real-time applications in automotive systems

##  II Overview

While our primary contribution is aggregation of the statistical method with a
Fully Connected Network (FCN), the former is also a novel method. This method
has been conceptualized by the first author and has been thoroughly adapted to
an industrial setting with the help of the second author, proving to produce
reliable results on its own. However, its full potential is realized when used
in conjunction with a deep learning model. The deep learning model makes
predictions and provides confidence scores. In cases where the confidence is
low, the statistical method is applied to aggregate the predictions and
produce a final, more accurate output.

The contributions of each author for the work presented in this paper are as
follows:

  * •

First author: The first author conceptualized the statistical approach. He is
also responsible for the creation of the artificial labeling, development of
the aggregation function for labels, and conduction of the statistical
method’s experiments.

  * •

Second author: The second author played a crucial role in transitioning the
theoretical method to an industrial application. Working with the baseline
method provided by the third author. Introduced the idea of SAAD and helped
showcase the results in this paper.

  * •

Third author: Provided the thoroughly tested baseline model. The author also
contributed by supplying and generating on demand the necessary data for
extensive experiments.

##  III Related Works

The field of anomaly detection has garnered significant attention due to the
increasing complexity and connectivity of modern systems. Several approaches
have been proposed to address the challenges associated with identifying
anomalies in various domains, with our particular interest in this paper being
the combination of machine learning and statistical methods.

Intrusion Detection Systems (IDS) for networks, particularly the Controller
Area Network (CAN) bus in automotive systems, have been extensively reviewed
and developed. Lokman et al. [19] provide a comprehensive overview of IDS
implementation in the CAN bus system, discussing various detection approaches,
deployment strategies, and technical challenges. They highlight the use of
hybrid methods that combine frequency-based, machine learning-based, and
statistical-based approaches for anomaly detection, emphasizing the importance
of a holistic approach to securing networks.

Statistical learning methods have also been applied to anomaly detection in
manufacturing environments, which share similarities with other systems in
terms of operational variability and data characteristics. Pittino et al. [20]
demonstrate the effectiveness of integrating statistical machine learning
methods to control charts in detecting anomalies in production machines. Their
hybrid approach, designed to function without continuous recalibration,
achieves high recall and robustness, showcasing the potential for similar
methods in various applications.

Deep learning techniques, particularly those involving Long Short-Term Memory
(LSTM) networks, have shown promise in handling the complex and unpredictable
nature of time-series data from multiple sensors. Malhotra et al. [21] propose
an LSTM-based Encoder-Decoder scheme for anomaly detection, which learns to
reconstruct normal time-series behavior and uses reconstruction error to
identify anomalies. This method, while primarily focused on machine learning,
can benefit from statistical preprocessing to enhance its effectiveness across
various types of time-series data, making it relevant for numerous
applications.

The combination of statistical methods with machine learning has been explored
in various fields, providing valuable insights for their application in
anomaly detection. Holloway and Mengersen [22] review the use of statistical
machine learning methods in remote sensing, highlighting their relevance for
sustainable development goals. Similarly, Amin et al. [23] demonstrate the
integration of statistical methods with machine learning for brain tumor
detection, achieving high accuracy through a combination of denoising
techniques and feature fusion.

ArunKumar et al. [24] utilize Auto-Regressive Integrated Moving Average
(ARIMA) and Seasonal Auto-Regressive Integrated Moving Average (SARIMA) models
to forecast COVID-19 cases, underscoring the importance of combining
statistical modeling with machine learning for accurate predictions.
Banitalebi Dehkordi et al. [25] apply a hybrid approach to detect DDoS attacks
in software-defined networks (SDN), integrating entropy-based methods with
machine learning classifiers to enhance detection accuracy.

In the context of cryptocurrency price prediction, Khedr et al. [26] review
traditional statistical and machine learning techniques, demonstrating the
effectiveness of machine learning in handling the dynamic and non-seasonal
nature of cryptocurrency data. Their work highlights the potential for similar
hybrid approaches in anomaly detection.

A structured approach to anomaly detection for in-vehicle networks is
presented by Muter et al. [27]. Their work focuses on a systematic methodology
for identifying anomalies in vehicular communications, integrating statistical
analysis with machine learning-based detection mechanisms to improve accuracy
and robustness.

Research Gap and Our Solution. Despite the significant advancements in anomaly
detection, there remains a gap in effectively integrating statistical methods
with machine learning models to achieve high accuracy and robustness,
particularly in real-time applications. Many existing methods either rely
heavily on machine learning techniques, which can be data-intensive and
computationally expensive, or on statistical methods that may not fully
capture the complexity of the data.

Our response to that is SAAD, which effectively demonstrates that standalone
models, which perform similarly when used individually, benefit greatly from
combining their results. By integrating statistical techniques with Fully
Connected Networks (FCNs), SAAD enhances the detection of anomalies by
leveraging the strengths of both approaches. Statistical preprocessing
improves data quality, while dropout layers in FCNs prevent overfitting,
leading to a more robust and accurate anomaly detection system.

##  IV Solution methodology

In this section, we detail the comprehensive methodology employed in the
development and validation of the SAAD framework. The methodology encompasses
three main stages: (1) the artificial labeling procedure, (2) the development
of the deep learning model, and (3) the aggregation of results from both
methods to achieve enhanced anomaly detection performance. For a visual
representation of the architecture of SAAD, we provide Figure 1.

###  IV-A Artificial labeling procedure

This stage primarily involves assigning artificial labels to dataset
instances. This approach is rooted in the fact that most anomaly detection
problems deal with unlabeled data. Typically, anomaly detection is an
unsupervised learning task due to the difficulty in defining anomalies and the
labor-intensive process of manually labeling data, which requires domain
expert knowledge. Starting with the initial dataset, but excluding the target
column, statistical methods are used to assign artificial or synthetic labels
to the data. This method offers a different perspective compared to the
supervised usage in deep learning, making their combination complementary. The
term “artificial” indicates that the labels generated in this stage are not
necessarily the ground-truth labels.

Anomalous measurements are identified for each feature of the dataset by
assessing the occurrence probabilities of the values. The concept originates
from the fact that anomalies are often not precisely defined; rather, they are
described as instances that exhibit the common behavior [28]. Following this
idea, in statistical terms, a value could be considered as being an anomaly if
its occurrence probability is lower than a given threshold value. Since it is
hard to work directly on high-dimensional data, each feature is examined
separately. The occurrence probability of a given value is determined using
the probability distribution function (PDF). Determining the exact probability
density function (PDF) is often challenging, making it more practical to
estimate the PDF using a histogram-based method. A histogram is created by
splitting the initial collection of numerical values into groups called bins.
Each bin has a defined lower and upper bound, which are used to determine
where each value should be placed. The number and widths of the bins could be
selected using different methods, but a proper width is determined using the
Freedman–Diaconis rule [29]. The number of bins is inferred from the usage of
bins with the same width.

The number of values from a given bin is called bin count. The relative
frequency of a bin is calculated by dividing the bin count by the total count
of all bins. Consider a bin being represented as
$b_{i}=\\{(a_{i-1},a_{i}),\>c_{i},\>m_{i}\\},\>i\geq 1$: $a_{i-1}$ and $a_{i}$
are the bin’s lower and upper bounds, $c_{i}$ is the bin count and $m_{i}$ is
the relative frequency. Bins with high counts contain significantly more
values, indicating that these values are more common within the initial
collection. Consequently, the values in high-count bins have higher occurrence
probabilities compared to those in low-count bins. A value that belongs to a
bin with a low count could be considered as being an anomaly because it is not
so common, i.e., it has a low number of occurrences. A bin $b_{i}$ is called
anomalous if its relative frequency value is lower than a given threshold
value $t$, i.e. $m_{i}<t$. An exemplification of non-anomalous and anomalous
bins on the histogram is presented in fig. 2. Further, a value
$x\in\mathbb{R}$ is considered as being an anomaly if it belongs to an
anomalous bin.

Figure 2: Anomalous bins exemplification on histogram

The selection of anomalous bins is presented in Algorithm 1. $B$ is a
collection of collections: each collection $B_{j}$ contains the histogram bins
for the values of feature $j$. The collection $A$ is used to retain the
anomalous bins for each feature and is structured in the same way as $B$.

Algorithm 1 Anomalous bins selection for all features

0: dataset $D\subseteq\mathbb{R}^{m\times n}$, $t\in[0,1]$, bin width $k$
(optional, use ’auto’ for Freedman–Diaconis rule)

0: $B=\\{B_{0},B_{1},...,B_{n-1}\\}$, $A=\\{A_{0},A_{1},...,A_{n-1}\\}$

$A$,
$B\leftarrow[{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}x}],[{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}x}]$
{empty lists}

$j\leftarrow 0$

for $j\leftarrow 0$ to $n-1$ do {for each feature}

$D^{j}\leftarrow$ values of feature j from all instances

$B_{j}\leftarrow get\\_bins(D^{j},\,k)$

add $B_{j}$ to $B$

$i\leftarrow 0$

$A_{j}\leftarrow[{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}x}]$
{empty list for anomalous bins}

for $i\leftarrow 0$ to $|B_{j}|-1$ do {for each bin}

$b_{i}\leftarrow B_{j,i}$ {get current bin}

if $m_{i}<t$ then {compare relative freq. with threshold}

add $b_{i}$ to $A_{j}$

end if

end for

add $A_{j}$ to $A$

end for

return $B$, $A$

After the selection of anomalous bins for each feature, an instance-level
anomaly is defined. An instance is considered an anomaly (label value 1) if it
contains at least $k$ anomalous values (features). The procedure is presented
in Algorithm 2.

Algorithm 2 Artificial labeling procedure

0: instance $d\in\mathbb{R}^{1\times n}$, $B=\\{B_{0},B_{1},...,B_{n-1}\\}$,
$A=\\{A_{0},A_{1},...,A_{n-1}\\}$, $k<n$

0: 1 if $d$ is an (artificial) anomaly, 0 otherwise

$count\leftarrow 0$

for $j\leftarrow 0$ to $n-1$ do

$b\leftarrow\texttt{bin\\_of\\_value}(d_{j},B_{j})$

if $b$ in $A_{j}$ then

$count\leftarrow count+1$

end if

end for

if $count\geq k$ then

return $1$

end if

return $0$

###  IV-B Deep Learning Model

The data preparation process is crucial for providing the model with clean and
well-structured input data. Initially, data is loaded from CSV files located
in a specified directory, generating file paths and reading the CSV files.
During loading, the data undergoes cleaning to handle missing or undefined
values, ensuring the dataset is free from any irregularities that could affect
model training. Relevant features, or signals, are then selected from the
dataset based on the specified configuration, focusing on the necessary
columns for model training. Categorical features have been encoded and
subsequently, the data is labeled based on the presence of anomalies, with a
binary label (0 or 1) assigned to each data point to indicate whether it is an
anomaly.

After labeling, the data is scaled using a standard scaler to normalize the
feature values. This normalization step helps in speeding up the convergence
during training by ensuring that all features contribute equally to the
learning process.

The dense network architecture comprises several fully connected (dense)
layers with dropout layers between them to prevent overfitting. The input
layer uses the selected relevant features and feeds them to the hidden layers,
which consist of multiple dense layers with ReLU activation functions and L2
regularization to mitigate overfitting. The output layer is a single neuron
with a sigmoid activation function to produce a binary output indicating the
presence or absence of an anomaly.

The model is compiled using the Adam optimizer, a binary cross entropy loss
function, and various evaluation metrics including binary accuracy, precision,
recall, and custom F-beta scores.

During inference, the model processes new input data in the same manner,
outputting a prediction along with a confidence score. The prediction is
binary, indicating whether the input data point is an anomaly, while the
confidence score reflects the model’s certainty about this prediction. The
confidence score is particularly important in applying SAAD, as the
aggregation of results from the deep learning model and the statistical model
is based on the confidence from the former.

###  IV-C Aggregated solution

The results provided by the artificial labeling procedure and the deep
learning model are aggregated to create a new set of labels. These final
results are compared with the ground-truth labels to compute the performance
metrics for the proposed approach. For an instance $d\in D$, consider
STAT_LABEL as the label (prediction result) provided by the artificial
labeling procedure and DL_LABEL as the label provided by the deep learning
model for the evaluation of $d$. For an instance $d$ evaluated by the deep
learning model, the outcome value DL_LABEL (either 1 or 0) depends on the
value $v$, computed at the end of the evaluation using the sigmoid function.
The value $v$ is further used by the deep learning model to decide the outcome
for the evaluated instance $d$: if, $v>0.5$ then, the model classifies $d$ as
an anomaly (label value 1), otherwise as a non-anomalous instance (label value
0). The aggregation function decides the final label value using the values of
STAT_LABEL, DL_LABEL, and $v$. The decision cases are described below:

STAT_LABEL | DL_LABEL | Condition | Result
---|---|---|---
0 | 0 | - | 0
1 | 1 | - | 1
0 | 1 | $v>a$ | 1
0 | 1 | $v\leq a$ | 0
1 | 0 | $v<b$ | 1
1 | 0 | $v\geq b$ | 0
TABLE I: Decision rules for label aggregation TABLE II: Accuracy Scores Obtained on Aggregated Results (in %) $b$ \$a$ | 0.51 | 0.52 | 0.53 | 0.54 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95 | 0.96 | 0.97 | 0.98 | 0.99 | 1
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
0.00 | 54.0 | 54.3 | 54.6 | 54.8 | 54.9 | 56.2 | 59.1 | 59.9 | 61.1 | 62.2 | 63.4 | 65.4 | 68.3 | 73.0 | 77.1 | 77.1 | 77.1 | 77.1
0.01 | 61.8 | 62.0 | 62.3 | 62.5 | 62.7 | 63.9 | 66.9 | 67.6 | 68.8 | 69.9 | 71.2 | 73.2 | 76.0 | 80.7 | 84.8 | 84.8 | 84.8 | 84.8
0.02 | 62.1 | 62.4 | 62.7 | 62.9 | 63.1 | 64.3 | 67.2 | 68.0 | 69.2 | 70.3 | 71.5 | 73.5 | 76.4 | 81.1 | 85.2 | 85.2 | 85.2 | 85.2
0.03 | 62.4 | 62.7 | 62.9 | 63.1 | 63.3 | 64.6 | 67.5 | 68.3 | 69.4 | 70.6 | 71.8 | 73.8 | 76.6 | 81.3 | 85.4 | 85.4 | 85.4 | 85.4
0.04 | 62.6 | 62.8 | 63.1 | 63.3 | 63.5 | 64.7 | 67.6 | 68.4 | 69.6 | 70.7 | 71.9 | 74.0 | 76.8 | 81.5 | 85.6 | 85.6 | 85.6 | 85.6
0.05 | 62.7 | 63.0 | 63.2 | 63.4 | 63.6 | 64.8 | 67.7 | 68.5 | 69.7 | 70.8 | 72.1 | 74.1 | 76.9 | 81.6 | 85.7 | 85.7 | 85.7 | 85.7
0.10 | 63.8 | 64.0 | 64.3 | 64.5 | 64.7 | 65.9 | 68.8 | 69.6 | 70.7 | 71.1 | 72.4 | 74.4 | 77.2 | 81.9 | 86.0 | 86.0 | 86.0 | 86.0
0.15 | 64.4 | 64.7 | 64.9 | 65.1 | 65.3 | 66.5 | 69.4 | 70.2 | 71.3 | 71.4 | 72.6 | 74.6 | 77.4 | 82.1 | 86.2 | 86.2 | 86.2 | 86.2
0.20 | 64.8 | 65.0 | 65.3 | 65.5 | 65.7 | 66.9 | 69.8 | 70.6 | 71.7 | 71.6 | 72.8 | 74.8 | 77.6 | 82.3 | 86.4 | 86.4 | 86.4 | 86.4
0.25 | 65.1 | 65.4 | 65.6 | 65.8 | 66.0 | 67.2 | 70.1 | 70.9 | 72.0 | 71.7 | 72.9 | 75.0 | 77.8 | 82.5 | 86.6 | 86.6 | 86.6 | 86.6
0.30 | 65.4 | 65.7 | 65.9 | 66.1 | 66.3 | 67.5 | 70.4 | 71.2 | 72.3 | 71.8 | 73.1 | 75.1 | 77.9 | 82.6 | 86.7 | 86.7 | 86.7 | 86.7
0.35 | 65.8 | 66.0 | 66.3 | 66.5 | 66.7 | 67.9 | 70.8 | 71.6 | 72.7 | 72.0 | 73.3 | 75.3 | 78.2 | 82.9 | 86.9 | 86.9 | 86.9 | 86.9
0.40 | 66.2 | 66.4 | 66.7 | 66.9 | 67.1 | 68.3 | 71.2 | 72.0 | 73.1 | 72.1 | 73.3 | 75.3 | 78.2 | 82.9 | 86.9 | 86.9 | 86.9 | 86.9
0.45 | 66.6 | 66.8 | 67.1 | 67.3 | 67.5 | 68.7 | 71.6 | 72.4 | 73.5 | 72.2 | 73.4 | 75.4 | 78.3 | 83.0 | 87.1 | 87.1 | 87.1 | 87.1
0.46 | 66.9 | 67.1 | 67.4 | 67.6 | 67.8 | 69.0 | 71.9 | 72.7 | 73.8 | 72.6 | 73.8 | 75.8 | 78.7 | 83.4 | 87.5 | 87.5 | 87.5 | 87.5
0.47 | 67.1 | 67.4 | 67.6 | 67.8 | 68.0 | 69.2 | 72.1 | 72.9 | 74.0 | 72.8 | 74.0 | 76.0 | 78.9 | 83.6 | 87.7 | 87.7 | 87.7 | 87.7
0.48 | 67.3 | 67.6 | 67.8 | 68.0 | 68.2 | 69.4 | 72.3 | 73.1 | 74.2 | 73.0 | 74.2 | 76.2 | 79.1 | 83.8 | 87.9 | 87.9 | 87.9 | 87.9
0.49 | 67.6 | 67.8 | 68.1 | 68.3 | 68.5 | 69.7 | 72.6 | 73.4 | 74.5 | 73.2 | 74.4 | 76.4 | 79.3 | 84.0 | 88.1 | 88.1 | 88.1 | 88.1
0.50 | 67.8 | 68.0 | 68.3 | 68.5 | 68.7 | 69.9 | 72.8 | 73.6 | 74.7 | 73.4 | 74.6 | 76.6 | 79.5 | 84.2 | 88.3 | 88.3 | 88.3 | 88.3
TABLE III: F1 Scores Obtained on Aggregated Results $b$ \$a$ | 0.51 | 0.52 | 0.53 | 0.54 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95 | 0.96 | 0.97 | 0.98 | 0.99 | 1
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
0.00 | 0.52 | 0.53 | 0.53 | 0.54 | 0.54 | 0.56 | 0.60 | 0.61 | 0.63 | 0.65 | 0.66 | 0.69 | 0.72 | 0.78 | 0.82 | 0.82 | 0.82 | 0.82
0.01 | 0.63 | 0.64 | 0.64 | 0.64 | 0.65 | 0.66 | 0.70 | 0.71 | 0.72 | 0.74 | 0.75 | 0.78 | 0.81 | 0.85 | 0.89 | 0.89 | 0.89 | 0.89
0.02 | 0.64 | 0.64 | 0.65 | 0.65 | 0.65 | 0.67 | 0.70 | 0.71 | 0.73 | 0.74 | 0.76 | 0.78 | 0.81 | 0.85 | 0.89 | 0.89 | 0.89 | 0.89
0.03 | 0.64 | 0.65 | 0.65 | 0.65 | 0.66 | 0.67 | 0.71 | 0.72 | 0.73 | 0.75 | 0.76 | 0.78 | 0.81 | 0.86 | 0.89 | 0.89 | 0.89 | 0.89
0.04 | 0.65 | 0.65 | 0.65 | 0.66 | 0.66 | 0.67 | 0.71 | 0.72 | 0.73 | 0.75 | 0.76 | 0.78 | 0.81 | 0.86 | 0.89 | 0.89 | 0.89 | 0.89
0.05 | 0.65 | 0.65 | 0.65 | 0.66 | 0.66 | 0.68 | 0.71 | 0.72 | 0.74 | 0.75 | 0.76 | 0.79 | 0.81 | 0.86 | 0.89 | 0.89 | 0.89 | 0.89
0.10 | 0.66 | 0.66 | 0.67 | 0.67 | 0.67 | 0.69 | 0.72 | 0.73 | 0.75 | 0.75 | 0.77 | 0.79 | 0.82 | 0.86 | 0.90 | 0.90 | 0.90 | 0.90
0.15 | 0.67 | 0.67 | 0.67 | 0.68 | 0.68 | 0.69 | 0.73 | 0.74 | 0.75 | 0.75 | 0.77 | 0.79 | 0.82 | 0.86 | 0.90 | 0.90 | 0.90 | 0.90
0.20 | 0.67 | 0.67 | 0.68 | 0.68 | 0.68 | 0.70 | 0.73 | 0.74 | 0.76 | 0.76 | 0.77 | 0.79 | 0.82 | 0.87 | 0.90 | 0.90 | 0.90 | 0.90
0.25 | 0.68 | 0.68 | 0.68 | 0.68 | 0.69 | 0.70 | 0.74 | 0.75 | 0.76 | 0.76 | 0.77 | 0.79 | 0.82 | 0.87 | 0.90 | 0.90 | 0.90 | 0.90
0.30 | 0.68 | 0.68 | 0.68 | 0.69 | 0.69 | 0.71 | 0.74 | 0.75 | 0.77 | 0.76 | 0.77 | 0.80 | 0.82 | 0.87 | 0.90 | 0.90 | 0.90 | 0.90
0.35 | 0.68 | 0.68 | 0.69 | 0.69 | 0.69 | 0.71 | 0.74 | 0.76 | 0.77 | 0.76 | 0.78 | 0.80 | 0.83 | 0.87 | 0.91 | 0.91 | 0.91 | 0.91
0.40 | 0.68 | 0.69 | 0.69 | 0.69 | 0.70 | 0.71 | 0.75 | 0.76 | 0.77 | 0.76 | 0.78 | 0.80 | 0.83 | 0.87 | 0.91 | 0.91 | 0.91 | 0.91
0.45 | 0.69 | 0.69 | 0.69 | 0.70 | 0.70 | 0.72 | 0.75 | 0.76 | 0.78 | 0.76 | 0.78 | 0.80 | 0.83 | 0.87 | 0.91 | 0.91 | 0.91 | 0.91
0.46 | 0.69 | 0.69 | 0.70 | 0.70 | 0.70 | 0.72 | 0.75 | 0.76 | 0.78 | 0.77 | 0.78 | 0.80 | 0.83 | 0.87 | 0.91 | 0.91 | 0.91 | 0.91
0.47 | 0.69 | 0.70 | 0.70 | 0.70 | 0.70 | 0.72 | 0.76 | 0.77 | 0.78 | 0.77 | 0.78 | 0.80 | 0.83 | 0.88 | 0.91 | 0.91 | 0.91 | 0.91
0.48 | 0.70 | 0.70 | 0.70 | 0.70 | 0.71 | 0.72 | 0.76 | 0.77 | 0.78 | 0.77 | 0.79 | 0.81 | 0.84 | 0.88 | 0.92 | 0.92 | 0.92 | 0.92
0.49 | 0.70 | 0.70 | 0.71 | 0.71 | 0.71 | 0.73 | 0.76 | 0.77 | 0.79 | 0.77 | 0.79 | 0.81 | 0.84 | 0.88 | 0.92 | 0.92 | 0.92 | 0.92
0.50 | 0.70 | 0.71 | 0.71 | 0.71 | 0.71 | 0.73 | 0.77 | 0.78 | 0.79 | 0.78 | 0.79 | 0.81 | 0.84 | 0.89 | 0.92 | 0.92 | 0.92 | 0.92

##  V Experiment

Figure 3: Heatmap of Accuracy Scores Obtained on Aggregated Results. The
heatmap visualizes the accuracy percentages across different thresholds $a$
and $b$. The color intensity represents the accuracy, with warmer colors
indicating higher accuracy scores. The pattern shows that accuracy generally
increases as both thresholds $a$ and $b$ increase, with the highest accuracy
observed in the lower right corner of the heatmap, where both thresholds are
at their maximum values. This indicates a positive correlation between higher
threshold values and improved accuracy.

###  V-A Data Acquisition

The data acquisition process for this study focuses on measuring and
collecting data from the Hardware-in-the-Loop (HIL) environment, specifically
tailored for underbrake scenarios. This is also presented in the Figure 1. The
workbench setup includes several components to simulate and capture the
necessary data:

Simulation Workbench. The simulation workbench connects a real-time simulation
environment to the physical Electronic Control Unit (ECU) of a vehicle. The
simulation environment encompasses detailed models of various vehicle systems
such as the engine, transmission, suspension, and steering. These models
interact with sensor readings, control signals, and actuator commands to mimic
real-world driving conditions.

The HIL setup provides inputs to the ECU, which then processes these inputs
using its control algorithms. The vehicle’s response, including feedback from
sensors such as wheel speed, brake pressure, and brake pad temperature is sent
back to the simulation environment. A benefit of this feedback loop is that it
allows engineers to test and validate the performance of the vehicle’s braking
system and identify any unwanted behavior in the ECU’s algorithms without the
need for physical testing.

Data Acquisition System (DAQ). The DAQ system is a critical component of the
workbench, designed to measure and collect high-fidelity data from the HIL
environment. The DAQ system comprises the following elements:

  * •

Sensors: Devices that measure various parameters during the simulation,
including temperature, position, voltage, and other relevant metrics.

  * •

Data Acquisition Hardware: This hardware interfaces the sensors with the
computer and typically includes an analog-to-digital converter to transform
the analog signals from the sensors into digital signals that can be processed
by the computer.

  * •

Software: The DAQ software is used to configure the hardware, process the
collected data, and display the results. It includes drivers necessary for
data acquisition and tools for data analysis.

  * •

Computer System: A robust computer system is utilized to run the DAQ software
and store the collected data for subsequent analysis.

Data Collection Process. The raw data is represented by measurements obtained
from the ECU within the HIL environment. Each measurement consists of
thousands of signals sampled at 10ms intervals, with each signal potentially
representing a unique feature.

To gather these measurements, various vehicle maneuvers were simulated to
reproduce underbrake scenarios. Key variables for these maneuvers included
vehicle speed (ranging from 0 to 150 km/h) and brake pedal pressure (ranging
from 5% to 100%). These maneuvers ensured a diverse and representative
dataset, which is essential for training robust machine learning models.

Experimental Setup. For this study, a total of 150 vehicle maneuvers were
simulated for the training set to capture a wide range of underbrake
scenarios. Additionally, distinct maneuvers were created for the testing set,
featuring multiple underbrake activations and varying speeds. This approach
helps to evaluate the performance of the anomaly detection algorithms and
prevents overfitting.

By carefully designing the data acquisition and preprocessing steps, we ensure
that the machine learning models trained on this data can make accurate
predictions and decisions, ultimately improving the reliability of automotive
systems.

###  V-B Experiment methodology

As proof of concept, the proposed method was used to detect anomalies on the
data collected from a HIL environment: as mentioned in previous section, the
used data represent measurements taken from different usages of braking
function of vehicles on different situations; the used datasets contain 18
features and the target class.

Models parameters configuration. In the case of artificial labeling procedure,
the used parameters are: number of minimum anomalous values per instance $k=3$
and anomalous bin relative frequency threshold $t=0.05$. The artificial
labeling method works on an unsupervised way: it assigns labels just by usages
of statistical methods, without using ground truth labels. The used parameters
values were chosen such that, the percentage of artificial anomalous labels to
be closer to the percentage of ground-truth anomalous labels.

For the deep learning model, to further enhance the model’s performance,
hyperparameter optimization is applied using Bayesian optimization with Keras
Tuner. The optimizer defines a search space for several key hyperparameters,
including the number of units in each dense layer, learning rate, dropout
rates, and regularization parameters.

The optimization process iteratively evaluates different combinations of these
hyperparameters by training the model on the training data and validating it
on a separate validation set. Early stopping is employed to prevent
overfitting, halting training when the model’s performance on the validation
set ceases to improve. The best-performing model configuration is then
selected based on the evaluation metrics.

###  V-C Experiments results

The artificial labeling method alone provided an accuracy score of 72.1% on
the testing data. The deep learning model alone provided on accuracy score of
71.5% on the testing data. The labels provided by each method were aggregated
according to the presented strategy. The new labels were compared with ground-
truth labels and accuracy and F1 score metrics were used to measure the
performances: the results are presented in the table II and III.

As can be observed from results from the tables II and III, the proposed
approach that combines artificial labeling method with deep learning model
provide higher scores than each method alone. Both parameters $a$ and $b$ from
aggregation function have an important role on the quality of the results.
Higher $a$ and $b$ allow the used approach to provide better results: labels
obtained with artificial labeling method together with higher $a$ and $b$
values enforce the aggregation function to select high confidence results
provided by deep learning model. The experimental results for different pairs
of $a$, $b$ values can also be seen in Figure 4.

Figure 4: Evolution of accuracy at different threshold values.

###  V-D Results interpretation

As presented in the table I, the aggregation approach provided higher
performances, in comparison with both models alone. As an initial explanation
for the linear / uniform increase in the accuracy and F1 score values, this
trend is likely related to the use of small steps in the values used in the
tuning for aggregation parameters $a$ and $b$. As explanations related to the
higher performances provided by the aggregation strategy, firstly for the case
when both models provide the same label, the aggregated label is well
understood.

When the deep learning model classifies an instance as an anomaly (label value
1) and the statistical approach classifies it as a normal instance (label
value 0), we rely on the deep learning model’s result only if it provides a
high enough confidence for the given prediction. For example, if $v=0.55$ this
value is close to 0.5, which is the upper bound for what is considered normal.
Thus, we use $v>a$ to ensure that the prediction quality is sufficient to
consider the evaluated instance as an anomaly. By increasing the value of $a$,
we request a higher confidence $v$ (expect progressively better result),
effectively raising the threshold for accepting its predictions. This higher
expectation leads to a decrease in the number of instances labeled as
anomalies (value 1) provided by the deep learning model. If $a=1$, the
condition $v>1$ holds nowhere, meaning all instances will be labeled as normal
(value 0), relying solely on the statistical approach’s labels.

When the deep learning model classifies an instance as non-anomalous (label
value 0) and the statistical approach classifies it as an anomaly (label value
1), we rely on the deep learning model if its confidence value $v$ is low
enough. A $v$ value lower than 0.5 but still close to it indicates a low
confidence level for the evaluated instance to be non-anomalous, so we enforce
prediction quality using a constraint rule based on parameter $b$: $v<b$. In
the case of the deep learning model, an instance $d$ is considered as being
normal if $v\in[0,0.5]$; thus, if $b=0.5$, the condition $v<0.5$ always holds,
thus all such cases receive the label 0 (not anomaly), according to the value
provided by the deep learning model.

To reiterate, for $a=1$ and $b=0.5$, when the two models provide distinct
labels for an evaluated instance, related to the outcomes provided by the
aggregation method, all anomaly labels (value 1) are coming from the
statistical approach, while all non-anomaly labels (value 0) come from the
deep learning model. The best performance is provided by the aggregation
function using the parameters $a=1$ and $b=0.5$ indicating that the
statistical approach performs better at predicting anomalies, whereas the deep
learning model excels at predicting non-anomalous instances. An explanation
for previous observation relates to the nature of the used models. The deep
learning model, due to the unbalanced nature of the dataset (with many more
normal instances compared to anomalies), may have a bias favoring the
prediction of normal instances over anomalies, leading to better performance
in predicting normal instances. Similarly, the statistical approach’s
superiority in predicting anomalies is likely due to its mode of operation: it
determines the anomalous or non-anomalous state of an evaluated instance by
examining feature values and checking for low-frequency values. This approach
uses several raw thresholds to decide what is an anomaly and what is normal.
The statistical approach’s outcome can heavily rely on just a few feature
values of the evaluated instance. For example, with $k=4$, means that if the
evaluated instance contains at least 4 anomalous values, the instance is
classified as an anomaly, even if the rest of the attributes have common
values (with respect to the PDF). Therefore, the statistical approach is
stricter regarding the values of an instance, making it easier to declare an
instance as anomalous: a sufficient number of anomalous values can override
the common (non-anomalous) values in the rest of the attributes.

##  VI Conclusion

In this paper, we introduced SAAD, a novel statistical anomaly detection
method with promising results for automotive systems, validated through
extensive testing using real-world data collected from a Hardware-in-the-Loop
(HIL) environment. Our main contribution is merging a novel statistical method
with deep learning, taking the best of both approaches and achieving higher
accuracy with the joint operation than with standalone models.

The results of our experiments demonstrate the efficacy of our proposed
method. The standalone statistical method achieved an accuracy of 72.1% on the
testing data, while the deep learning model alone provided an accuracy of
71.5%. When combined, the aggregated results yielded an accuracy of up to
88.3% and an F1 score of 0.921.

The integration of statistical techniques with machine learning models
highlights the potential for enhancing existing anomaly detection frameworks,
making them more reliable and effective in real-world applications.

Future work will focus on further refining the statistical method, exploring
adaptive threshold techniques for multivariate time series data, and expanding
the application scope to other critical automotive systems. This will help
improve the robustness and scalability of our approach, ensuring its
effectiveness across diverse and complex scenarios, while also extending its
applicability beyond automotive systems.

## Acknowledgment

We extend our gratitude to Prof. Darian Onchis for his invaluable guidance and
support during the planning stages of this work, and for shaping the direction
of our research. His initial works on hybrid models were crucial in
recognizing their potential and applicability in the industry.

## References

  * [1] S. A. Jebur, K. A. Hussein, H. K. Hoomod, L. Alzubaidi, and J. Santamaría, “Review on Deep Learning Approaches for Anomaly Event Detection in Video Surveillance,” _Electronics_ , vol. 12, no. 1, p. 29, Jan. 2023, number: 1 Publisher: Multidisciplinary Digital Publishing Institute. [Online]. Available: https://www.mdpi.com/2079-9292/12/1/29
  * [2] A. A. Alsulami, Q. Abu Al-Haija, A. Alqahtani, and R. Alsini, “Symmetrical Simulation Scheme for Anomaly Detection in Autonomous Vehicles Based on LSTM Model,” _Symmetry_ , vol. 14, no. 7, p. 1450, Jul. 2022, number: 7 Publisher: Multidisciplinary Digital Publishing Institute. [Online]. Available: https://www.mdpi.com/2073-8994/14/7/1450
  * [3] P. Dixit, P. Bhattacharya, S. Tanwar, and R. Gupta, “Anomaly detection in autonomous electric vehicles using ai techniques: A comprehensive survey,” _Expert Systems_ , vol. 39, no. 5, p. e12754, 2022.
  * [4] Q. He, “A machine learning-based anomaly detection framework for connected and autonomous vehicles cyber security,” Ph.D. dissertation, University of Nottingham, 2021.
  * [5] B. Kim, M. A. Alawami, E. Kim, S. Oh, J. Park, and H. Kim, “A Comparative Study of Time Series Anomaly Detection Models for Industrial Control Systems,” _Sensors_ , vol. 23, no. 3, p. 1310, Jan. 2023, number: 3 Publisher: Multidisciplinary Digital Publishing Institute. [Online]. Available: https://www.mdpi.com/1424-8220/23/3/1310
  * [6] G. Matieş and C. Fosalau, “Detection of abs events in electronic brake systems using machine learning algorithms,” in _2023 13th International Symposium on Advanced Topics in Electrical Engineering (ATEE)_. IEEE, 2023, pp. 1–6.
  * [7] A. Huet, J. M. Navarro, and D. Rossi, “Local Evaluation of Time Series Anomaly Detection Algorithms,” in _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , ser. KDD ’22. New York, NY, USA: Association for Computing Machinery, Aug. 2022, pp. 635–645. [Online]. Available: https://doi.org/10.1145/3534678.3539339
  * [8] N. Aslam, I. U. Khan, A. Alansari, M. Alrammah, A. Alghwairy, R. Alqahtani, R. Alqahtani, M. Almushikes, and M. A. Hashim, “Anomaly Detection Using Explainable Random Forest for the Prediction of Undesirable Events in Oil Wells,” _Applied Computational Intelligence and Soft Computing_ , vol. 2022, p. e1558381, Aug. 2022, publisher: Hindawi. [Online]. Available: https://www.hindawi.com/journals/acisc/2022/1558381/
  * [9] Y. Kim and H. K. Kim, “Anomaly Detection using Clustered Deep One-Class Classification,” in _2020 15th Asia Joint Conference on Information Security (AsiaJCIS)_ , Aug. 2020, pp. 151–157. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/9194140
  * [10] A. A. Cook, G. Mısırlı, and Z. Fan, “Anomaly Detection for IoT Time-Series Data: A Survey,” _IEEE Internet of Things Journal_ , vol. 7, no. 7, pp. 6481–6494, Jul. 2020, conference Name: IEEE Internet of Things Journal. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/8926446
  * [11] A. Blázquez-García, A. Conde, U. Mori, and J. A. Lozano, “A Review on Outlier/Anomaly Detection in Time Series Data,” _ACM Computing Surveys_ , vol. 54, no. 3, pp. 56:1–56:33, Apr. 2021. [Online]. Available: https://doi.org/10.1145/3444690
  * [12] H. D. M. Onchis, “A deep learning approach to condition monitoring of cantilever beams via time-frequency extended signatures,” _Computers in Industry_ , vol. 105, pp. 177–181, Feb. 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0166361518305669
  * [13] D. M. Onchis and G.-R. Gillich, “Stable and explainable deep learning damage prediction for prismatic cantilever steel beam,” _Computers in Industry_ , vol. 125, p. 103359, Feb. 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0166361520305935
  * [14] D. Onchis, C. Istin, and E. Hogea, “A Neuro-Symbolic Classifier with Optimized Satisfiability for Monitoring Security Alerts in Network Traffic,” _Applied Sciences_ , vol. 12, no. 22, p. 11502, Jan. 2022, number: 22 Publisher: Multidisciplinary Digital Publishing Institute. [Online]. Available: https://www.mdpi.com/2076-3417/12/22/11502
  * [15] D. M. Onchis, G.-R. Gillich, E. Hogea, and C. Tufisi, “Neuro-symbolic model for cantilever beams damage detection,” _Computers in Industry_ , vol. 151, p. 103991, Oct. 2023. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0166361523001410
  * [16] F. Costi, D. Onchis, E. Hogea, and C. Istin, “Predictive Modeling for Diabetes Using GraphLIME,” Mar. 2024, pages: 2024.03.14.24304281. [Online]. Available: https://www.medrxiv.org/content/10.1101/2024.03.14.24304281v1
  * [17] E. Hogea, A. Popescu, D. Onchis, and G. Petit, “FeTrIL++: Feature Translation for Exemplar-Free Class-Incremental Learning with Hill-Climbing,” Mar. 2024, arXiv:2403.07406 [cs]. [Online]. Available: http://arxiv.org/abs/2403.07406
  * [18] D. M. Onchis, C. Istin, and H. Eduard-Florin, “Advantages of a neuro-symbolic solution for monitoring IT infrastructures alerts,” in _2022 24th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC)_. Hagenberg / Linz, Austria: IEEE, Sep. 2022, pp. 189–194. [Online]. Available: https://ieeexplore.ieee.org/document/10130897/
  * [19] S.-F. Lokman, A. T. Othman, and M.-H. Abu-Bakar, “Intrusion detection system for automotive Controller Area Network (CAN) bus system: a review,” _J Wireless Com Network_ , vol. 2019, no. 1, p. 184, Dec. 2019. [Online]. Available: https://jwcn-eurasipjournals.springeropen.com/articles/10.1186/s13638-019-1484-3
  * [20] F. Pittino, M. Puggl, T. Moldaschl, and C. Hirschl, “Automatic Anomaly Detection on In-Production Manufacturing Machines Using Statistical Learning Methods,” _Sensors_ , vol. 20, no. 8, p. 2344, Jan. 2020, number: 8 Publisher: Multidisciplinary Digital Publishing Institute. [Online]. Available: https://www.mdpi.com/1424-8220/20/8/2344
  * [21] P. Malhotra, A. Ramakrishnan, G. Anand, L. Vig, P. Agarwal, and G. Shroff, “LSTM-based Encoder-Decoder for Multi-sensor Anomaly Detection,” Jul. 2016, arXiv:1607.00148 [cs, stat]. [Online]. Available: http://arxiv.org/abs/1607.00148
  * [22] J. Holloway and K. Mengersen, “Statistical Machine Learning Methods and Remote Sensing for Sustainable Development Goals: A Review,” _Remote Sensing_ , vol. 10, no. 9, p. 1365, Sep. 2018, number: 9 Publisher: Multidisciplinary Digital Publishing Institute. [Online]. Available: https://www.mdpi.com/2072-4292/10/9/1365
  * [23] J. Amin, M. Sharif, M. Raza, T. Saba, and M. A. Anjum, “Brain tumor detection using statistical and machine learning method,” _Computer Methods and Programs in Biomedicine_ , vol. 177, pp. 69–79, Aug. 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0169260718313786
  * [24] K. E. ArunKumar, D. V. Kalaga, C. M. Sai Kumar, G. Chilkoor, M. Kawaji, and T. M. Brenza, “Forecasting the dynamics of cumulative COVID-19 cases (confirmed, recovered and deaths) for top-16 countries using statistical machine learning models: Auto-Regressive Integrated Moving Average (ARIMA) and Seasonal Auto-Regressive Integrated Moving Average (SARIMA),” _Applied Soft Computing_ , vol. 103, p. 107161, May 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1568494621000843
  * [25] A. Banitalebi Dehkordi, M. Soltanaghaei, and F. Z. Boroujeni, “The DDoS attacks detection through machine learning and statistical methods in SDN,” _J Supercomput_ , vol. 77, no. 3, pp. 2383–2415, Mar. 2021. [Online]. Available: https://doi.org/10.1007/s11227-020-03323-w
  * [26] A. M. Khedr, I. Arif, P. R. P V, M. El-Bannany, S. M. Alhashmi, and M. Sreedharan, “Cryptocurrency price prediction using traditional statistical and machine-learning techniques: A survey,” _Intelligent Systems in Accounting, Finance and Management_ , vol. 28, no. 1, pp. 3–34, 2021, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/isaf.1488. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/isaf.1488
  * [27] M. Muter, A. Groll, and F. C. Freiling, “A structured approach to anomaly detection for in-vehicle networks,” in _2010 Sixth International Conference on Information Assurance and Security_. Atlanta, GA, USA: IEEE, Aug. 2010, pp. 92–98. [Online]. Available: http://ieeexplore.ieee.org/document/5604050/
  * [28] L. Ruff, J. R. Kauffmann, R. A. Vandermeulen, G. Montavon, W. Samek, M. Kloft, T. G. Dietterich, and K.-R. Muller, “A unifying review of deep and shallow anomaly detection,” _Proceedings of the IEEE_ , vol. 109, no. 5, p. 756–795, May 2021. [Online]. Available: http://dx.doi.org/10.1109/JPROC.2021.3052449
  * [29] D. Freedman and P. Diaconis, “On the histogram as a density estimator: L2 theory,” _Probability Theory and Related Fields_ , vol. 57, no. 4, pp. 453–476, Dec 1981. [Online]. Available: https://doi.org/10.1007/BF01025868


††thanks: Phys. Rev. D 104 (2021) 076004,
https://link.aps.org/doi/10.1103/PhysRevD.104.076004

# Reinterpretation of proton strangeness between ATLAS and CMS measurements

Alim Ruzi  Bo-Qiang Ma<EMAIL_ADDRESS>
###### Abstract

We reexamine the shapes of the strange quark parton distribution functions
(PDFs) of the proton by means of QCD analysis of HERA I+II $e^{\pm}p$ deep
inelastic scattering cross section measurement at DESY, and inclusive gauge
boson production and $W$ boson production associated with a charm quark from
the LHC at CERN. We find that there is an overall agreement on the strange
quark distributions obtained from CMS $W$ \+ charm and ATLAS $W/Z$ data at the
parton momentum fraction range $x\lesssim 10^{-2}$. Meanwhile, there is also a
strong tension between these data toward large $x$. We find that this tension
fades away if the ATLAS measurement of $W/Z$ production is analyzed together
with the ATLAS $W$ \+ charm data. The $W/Z$ and $W$ \+ charm data both from
ATLAS and CMS experiments agree that the proton strangeness is enhanced toward
small momentum fraction $x$ and is smoothly suppressed at large $x$.
Furthermore, a strong $x$ dependence of the strange-to-nonstrange parton ratio
$R_{s}(x,Q^{2})$ is observed.

proton, strangeness, parton distribution function, deep inelastic scattering
proton, proton collision, inclusive gauge boson production

##  I Introduction

It is well known that nucleons consist of pointlike particles denominated as
partons, revealed by lepton nucleon deep inelastic scattering (DIS) experiment
[1, 2] almost half a century ago. There have been tremendous efforts to
precisely determine the parton distribution functions (PDFs) of the proton
worldwide. Distributions of the valence-up quarks, valence-down quarks and
gluons are determined to a level that could give, by and large, a consistent
result on physical quantities in the standard model (SM) obtained from high-
energy collision experiments [3, 4, 5, 6, 7]. However, the strange content of
the proton, which is poorly known because of limited sources of strangeness-
sensitive data, is still ambiguous to our understandings of the partonic
structure of the proton. Strange quark PDFs of nucleons play a vital role for
a number of physics processes at the LHC, ranging from measurements of $W$
boson production in association with a charm jet [8] and determinations of
electroweak (EW) interaction parameters to the formation of strange hadrons
[9].

The size and shape of the strange quark PDFs have recently attracted a lot of
interest and ignited a hot debate on specialized strange studies. The strange
quark PDFs obtained from dimuon production in neutrino-nucleus DIS experiment
by the NuTeV and NOMAD Collaborations [10, 11, 12], are significantly
suppressed while the inclusive $W/Z$ production data by the ATLAS
Collaboration [13] give strongly enhanced distribution of strangeness. The
ratio of strange to nonstrange sea quark PDFs, $R_{s}$ =
$(s+\bar{s})/(\bar{u}+\bar{d})$, from neutrino-nucleus DIS experiment is
approximately 0.5 when evaluated at $x$ = 0.023 and energy scale $Q$ = 1.6
GeV, whereas the strange to antidown quark ratio,
$r_{s}=(s+\bar{s})/2\bar{d}$, evaluated at $x$ = 0.023 and $Q^{2}$ = 1.9
$\mathrm{GeV^{2}}$ from the ATLAS analysis of the inclusive gauge boson
production together with the combined HERA I+II $e^{\pm}p$ data [14], rises up
to unity.

In their QCD analysis of accumulated $W/Z$ data to extract ATLAS-epWZ16 [15],
the ATLAS Collaboration obtains even larger strange-to-nonstrange ratio,
$R_{s}$ = 1.131 [15], showing agreement with their older analysis [13].
Another study of the ATLAS and CMS $W/Z$ cross section measurements in Ref.
[16] reports that there is no tension between the LHC $W/Z$ and the HERA data
or between the ATLAS and CMS datasets, indicating that the LHC measurements of
inclusive $W/Z$ boson production support unsuppressed strangeness in the
proton at low $x$ at both low-and high-energy scales.

A recent global analysis of vector gauge boson production in association with
jets plus older $W/Z$ and HERA $e^{\pm}$ data by the ATLAS Collaboration [17]
reports somewhat different $R_{s}$ distribution from the one in ATLAS-epWZ16
PDF sets. The ratio $R_{s}$ in [17] keeps consistency with ATLAS-epWZ16 at
$x\lesssim 10^{-2}$ and declines faster at other $x$ value. But, the strange
quark in that analysis still stays comparable or even larger than the
nonstrange sea quarks. The very recent global PDFs, MHST20 PDFs [18],
extracted from almost all of the available measurements show an $R_{s}$ value
compatible with that of the ATLAS-epWZ16.

So it comes as a puzzle that both the ATLAS and CMS measurements of $W/Z$
boson production support an enhanced strangeness inside the proton, while the
fixed target lepton-nuclear collision prefers a rather suppressed strangeness.
Among all the standard model processes, the $Ws$ $\to$ $c$ interactions have
the utmost sensitivity to the strange (antistrange) quark PDFs inside the
proton. The QCD analysis of $W$ production in association with a charm
(anticharm) quark in the proton-proton (pp) collision at $\sqrt{s}$ = 7 TeV by
the ATLAS Collaboration shows the ratio of strange to down sea quark,
$r_{s}=(s+\bar{s})/2\bar{d}$, to be 0.96 [19], supporting the results obtained
from the analysis of $W/Z$ measurements. However, the CMS Collaboration, in
the analysis of the same process, $W$ \+ charm production plus $W$-lepton
charge asymmetry [20, 21, 22, 23], reports that they have not observed the
enhanced strangeness as the ATLAS Collaboration did. Instead, they observed an
$x$-dependent and suppressed ratio $R_{s}$. Yet, in another dedicated analysis
on the strange quark distribution by the Hessian error updating method [24,
25], it is reported that the CMS $W$ \+ charm quark production data [22]
enhance strangeness of CT14 PDFs families. The recently extracted CT18 PDFs
series from fixed target and up-to-date collider data by the CTEQ group [26]
demonstrate that the strange quark density is enhanced over the nominal fit
CT18NNLO when the inclusive $W/Z$ data are included. Another devoted study on
the strangeness of the proton [27], which includes both the lepton-nucleus DIS
cross section data from the NuTeV and NOMAD plus $W$ \+ charm data from the
ATLAS and CMS experiments, reports a slightly enhanced strange quark PDF
within uncertainties in terms of the ratio $R_{s}(x,Q^{2})$ with huge
uncertainties at large momentum fraction.

Given the important influence of the strange quark PDFs in SM physics inside
the LHC, and the contradictory interpretations of the strange-sensitive data
from the ATLAS and CMS measurements, we regard reinvestigation of the exact
shapes of the strange quark PDFs of the nucleon by a global analysis of
strangeness-sensitive measurements, within an accurate theoretical and
methodological framework, as necessary to clarify the controversial
interpretations of the strangeness-related quantities, such as $xs(x,Q^{2})$
and the ratio $R_{s}$. This paper is intended to fulfill this purpose. We
present a fairly detailed PDFs extraction with the general-purpose PDFs
fitting program XFITTER [28] from the following data: HERA I+II $e^{\pm}p$
[14], ATLAS combined inclusive $W/Z$ production cross section [13, 15], and
$W$ production in association with a charm quark [19]; CMS $W$-lepton charge
asymmetry [20, 21] and $W$ production in association with charm quark [22,
23], applying variants of different parametrization forms first to clarify
whether the same physics processes in ATLAS and CMS experiments give
compatible or incompatible distributions of the strange quark and, second, to
pin down its exact shape.

##  II Input data and corresponding theoretical predictions

###  II.1 Data description

The process of $e^{\pm}p$ DIS has been the central dataset to many PDFs
working groups in their global analyses and provides the most constraining
power to quark PDFs when parametrized as a composite functional form. The
gauge boson production data at the LHC also provide significant constraints on
the flavor decomposition where the HERA data fail. The $W$ \+ charm production
offers the best constraints on the strange quark PDFs and gluon distributions
through the subprocesses: $g+s$ $\to$ $W^{-}$ \+ $c$ and $\bar{s}+g$ $\to$
$W^{+}$ \+ $\bar{c}$. The $W$-lepton charge asymmetry data from CMS, available
for $7~{}$ and $8~{}\mathrm{TeV}$, can constrain the valence quarks very well
and have an indirect impact on the strange content.

We do not include neutrino-induced dimuon data in our analysis. The dimuon
data were obtained from $\nu_{\mu}(\bar{\nu}_{\mu})A$ fixed target DIS process
where the target is iron, whose internal structure is more complex than the
free nucleon. In this process, however, the nuclear European Muon
Collaboration (EMC) effect should be taken into account in addition to some
nuclear corrections; see Ref. [29] and references therein. Considering the
non-negligible nuclear corrections in the lepton-nuclear collision experiment,
we believe that pure-proton-related data, which are free of the EMC effect,
are significantly useful for our analysis. Besides that, the strange quark
distribution extracted from dimuon production in neutrino scattering differs
significantly from the one extracted from LHC data, especially the ATLAS $W/Z$
production process. The latter prefers a rather higher $R_{s}$ value and the
former prefers a lower value. One more disadvantage of dimuon data is that
their constraining power vanishes at $x\lesssim 0.01$ due to shadowing and
anti-shadowing effects in neutrino-nucleus DIS process [30].

The description of data is given below:

  1. (i)

The HERA I+II $e^{\pm}p$ combined data [14] include measurements acquired from
run 1 during 1992 to 2000, and run 2 during 2002 to 2007 in $e^{\pm}p$
collision at different electron/positron and proton beam energies with
approximately the total integrated luminosities of 500 pb-1. The datasets are
divided between $e^{\pm}p$ charged current (CC) and neutral current (NC)
scattering according to the intermediate particle being $W^{\pm}$ or
$Z/\gamma^{*}$. The $e^{\pm}p$ CC and $e^{-}p$ NC interaction scattering cross
sections are measured at $\sqrt{s}$ = 318 GeV. The rest of the datasets are
made up of $e^{+}p$ NC scattering cross sections with positron beam energy of
$\mathrm{E_{e}}$ = 27.5 GeV and proton beam energies of $\mathrm{E_{p}}$ =
920, 820, 575 and 460 GeV, corresponding to total center-mass energies
$\sqrt{s}$ = 318, 300, 251, and 225 GeV. These datasets when combined give in
total 1307 data points on differential cross sections being a function of the
Bjorken variable $x$, approximately the light-cone momentum fraction of the
quark, and the transferred-energy squared $Q^{2}$. For the published datasets,
these two kinematic variables can cover $6\cdot 10^{-7}\leq x\leq 0.65$,
$0.045\leq Q^{2}\leq 50000$ $\mathrm{GeV^{2}}$ in the NC interaction
scattering and $1.3\cdot 10^{-2}\leq x\leq 0.40$, $200\leq Q^{2}\leq 50000$
$\mathrm{GeV^{2}}$ in the CC interaction scattering.

  2. (ii)

The ATLAS $W/Z$ data [15] include the cross section measurements of $W^{\pm}$
production in the leptonic decay channels: $W^{+}\to l^{+}\nu$, $W^{-}\to
l^{-}\bar{\nu}$ and $Z\to ll$ ($l=e,\mu$) production processes at $\sqrt{s}=$
7 TeV with an integrated luminosity of $4.6$ fb-1. The differential cross
sections for $W^{\pm}$ are given as the function of lepton pseudorapidity in
the range $|\eta_{l}|<$ 2.5 with 22 data points in total. The cross section of
$Z/\gamma^{*}$ production is measured as a function of the absolute dilepton
rapidity $|y_{ll}|$ in the central and forward region, for three intervals of
dilepton invariant mass or the intermediate $Z$ boson mass, $46<m_{ll}<66$
GeV, $66<m_{ll}<116$ GeV (forward and central rapidity), and $116<m_{ll}<150$
GeV (forward and central rapidity). The $Z/\gamma^{*}$ data include 39 data
points in total.

  3. (iii)

The data of $W$ \+ charm measurement come from ATLAS and CMS experiments. The
ATLAS measurements of $W^{\pm}$ production in association with an (anti)charm
quark are available for $\sqrt{s}$ = 7 TeV with total 22 data points [19],
while the CMS data of this process are available for 7 [22] and 13 TeV [23]
averaged for $W^{+}$ \+ $\bar{c}$ and $W^{-}$ \+ c with a total of five
points. Both ATLAS and CMS measurements of $W$ \+ charm production are given
as function of lepton pseudorapidity, $\eta^{\mu}$.

  4. (iv)

The $W$-lepton charge asymmetry, defined as $\mathcal{A}$ =
$\frac{\sigma(W^{+})-\sigma(W^{-})}{\sigma(W^{+})+\sigma(W^{-})}$, from CMS is
available for 7 [20] and 8 TeV [21].

###  II.2 Theoretical calculation

The corresponding predictions of $e^{\pm}p$ cross sections for HERA I+II
$e^{\pm}p$ data [14] are obtained by solving the DGLAP evolution equations
[31, 32, 33, 34] at next-to-next-to-leading-order (NNLO) in the
$\overline{\mathrm{MS}}$ scheme [35] through the QCDNUM program [36]
interfaced with XFITTER. The renormalization and factorization scales are
chosen to be $\mu_{r}^{2}$ = $\mu_{f}^{2}$ = $Q^{2}$. The initial
factorization scale is set to be $Q_{0}^{2}=1.9~{}\mathrm{GeV}^{2}$ just below
the charm quark mass $m_{c}=1.43~{}\mathrm{GeV}$ in this analysis. The evolved
PDFs at each scale point are convoluted with the coefficient functions of the
structure functions of the proton to obtain the corresponding differential
cross sections for $e^{\pm}p$ scattering. The contributions to the structure
functions from heavy quarks are calculated in the general-mass-variable-
flavor-number scheme [37, 38, 39] in the NC interactions. For the structure
functions in CC interactions, the zero-mass approximation is used because the
HERA CC data have $Q^{2}\gg M_{H}^{2}$, where $M_{H}$ is the heavy quark mass,
i.e., the masses of charm and bottom quarks. The calculation of the
differential cross section for the gauge boson production is done using the
well-known Monte Carlo program MCFM [40, 41] with the CT18NNLO PDF set [26].
The $W^{\pm}$ and $Z$ boson cross sections are available up to NNLO. Because
of the intensiveness of this calculation, a fast convolution technique is
applied. First of all, the partonic cross sections in $pp$ collision are
obtained in the form of grid files using MCFM with the help of another fast
parton convolution program APPLGRID[42] interfaced both with MCFM and XFITTER;
then, the corresponding cross sections are obtained convoluting the grid files
with the evolved PDFs during the fit. So far, APPLGRID (version-1.5.34) can
only calculate cross sections up to next-to-leading-order (NLO) accuracy. To
get the NNLO total and differential cross sections, a $K$-factor formalism is
applied. The $K$-factor is defined as

$K=\frac{\sigma^{\mathrm{NLO~{}EW}}_{\mathrm{NNLO~{}QCD}}}{\sigma^{\mathrm{LO~{}EW}}_{\mathrm{NLO~{}QCD}}}$
(1)

For both the CMS and ATLAS $W$ \+ charm data, only the NLO calculations are
used and this is obtained via MCFM. The theoretical prediction for $W$ \+
charm process is only available at NLO calculation, while that for the Drell-
Yan (DY) and $W^{\pm}$ production can be obtained up to NNLO. The measurements
for the $W/Z$ data [15] are given in high precision, so it is important to get
as accurate theoretical predictions as possible. During the fit, the order of
calculation is specified according to the physics process under consideration.
All the predictions are calculated in the respective fiducial phase space of
the experimental data. The $K$-factors are evaluated bin by bin with the same
PDFs in both the numerator and denominator.

##  III Analysis setup

We are curious about the ratio $R_{s}$ given by CMS analysis of $W$ \+ charm
production together with $W$-lepton charge asymmetry data [23] being much
different from the one reported by the ATLAS group [15]. The CMS analysis [23]
made a comparison of the ratio $R_{s}$ to the one obtained with ATLAS-epWZ16
PDFs [15] and ABMP16 PDFs [43]. $R_{s}$ from ATLAS is almost constant
throughout the entire $x$ except at high $x$ and the property of being greater
than 1 is interpreted as enhancement of strange quark PDFs relative to
nonstrange quark PDFs by the ATLAS group. On the contrary, $R_{s}$ of CMS is
fully $x$ dependent and declining faster as $x$ getting larger, interpreted as
suppressed strangeness by CMS. This is the main difference between ATLAS $W/Z$
and CMS $W$ \+ charm analyses in terms of the ratio $R_{s}$. Considered as
dominant data sources in constraining the light flavor quark and antiquark
PDFs, the ATLAS $W/Z$ and CMS $W$ \+ charm data should not be controversial as
far as the strange quark PDFs are concerned. Below, we are going to explore
both the analysis setup of the ATLAS [15] and CMS [23] analyses. The ATLAS-
epWZ16 is extracted at NNLO accuracy in perturbative QCD using inclusive HERA
I+II $e^{\pm}p$ [14] and ATLAS $W/Z$ data [15] through XFITTER framework. The
PDFs in ATLAS-epWZ16 are parametrized at a starting scale $Q^{2}_{0}$ = 1.9
$\mathrm{GeV}^{2}$ as

$\displaystyle xg(x)$
$\displaystyle=A_{g}x^{B_{g}}(1-x)^{C_{g}}-A^{\prime}_{g}x^{B^{\prime}_{g}}(1-x)^{C^{\prime}_{g}},$
(2a) $\displaystyle xu_{v}(x)$
$\displaystyle=A_{u_{v}}x^{B_{u_{v}}}(1-x)^{C_{u_{v}}}(1+E_{u_{v}}x^{2}),$
(2b) $\displaystyle xd_{v}(x)$
$\displaystyle=A_{d_{v}}x^{B_{d_{v}}}(1-x)^{C_{d_{v}}},$ (2c) $\displaystyle
x\bar{u}(x)$ $\displaystyle=A_{\bar{u}}x^{B_{\bar{u}}}(1-x)^{C_{\bar{u}}},$
(2d) $\displaystyle x\bar{d}(x)$
$\displaystyle=A_{\bar{d}}x^{B_{\bar{d}}}(1-x)^{C_{\bar{d}}},$ (2e)
$\displaystyle x\bar{s}(x)$
$\displaystyle=xs(x)=A_{\bar{s}}x^{B_{\bar{s}}}(1-x)^{C_{\bar{s}}}.$ (2f)

The relationship $A_{\bar{u}}$ = $A_{\bar{d}}$ and $B_{\bar{u}}$ =
$B_{\bar{d}}$ = $B_{\bar{s}}$ are set for the sea quark parameters based on
the assumption $\bar{u}$ = $\bar{d}$ as $x\to 0$, complying with the usual way
quoted by the HERA Collaboration [14].

On the contrary, the CMS Collaboration applies similar shape functions at
$Q^{2}_{0}$ = 1.9 $\mathrm{GeV^{2}}$ for their analysis of $W$ \+ charm quark
production together with the $W$-lepton charge asymmetry data to extract PDFs
at NLO using XFITTER. The applied parametrization form in the CMS analysis is:

$\displaystyle xg(x)$ $\displaystyle=A_{g}x^{B_{g}}(1-x)^{C_{g}},$ (3a)
$\displaystyle xu_{v}(x)$
$\displaystyle=A_{u_{v}}x^{B_{u_{v}}}(1-x)^{C_{u_{v}}}(1+E_{u_{v}}x^{2}),$
(3b) $\displaystyle xd_{v}(x)$
$\displaystyle=A_{d_{v}}x^{B_{d_{v}}}(1-x)^{C_{d_{v}}},$ (3c) $\displaystyle
x\bar{u}(x)$
$\displaystyle=A_{\bar{u}}x^{B_{\bar{u}}}(1-x)^{C_{\bar{u}}}(1+E_{\bar{u}}x^{2}),$
(3d) $\displaystyle x\bar{d}(x)$
$\displaystyle=A_{\bar{d}}x^{B_{\bar{d}}}(1-x)^{C_{\bar{d}}},$ (3e)
$\displaystyle x\bar{s}(x)$
$\displaystyle=xs(x)=A_{\bar{s}}x^{B_{\bar{s}}}(1-x)^{C_{\bar{s}}}.$ (3f)

In their parametrization form, the CMS Collaboration applies no extra
relationship between sea quark parameters except for the renormalization
parameters of $\bar{u}$ and $\bar{d}$, $A_{\bar{u}}=A_{\bar{d}}$. They put an
extra parameter that directly affects the strange quark PDFs: strangeness
fraction $f_{s}=\bar{s}/(\bar{d}+\bar{s})$, often used by the HERA group in
their PDFs analysis [14] as well. This strangeness fraction number is set as a
free parameter in the fit. Power parameters $B$ of all the light flavor sea
quarks, $\bar{u}$, $\bar{d}$, $\bar{s}$ are also set free from each other.

To know whether there is a tension between ATLAS $W/Z$ and CMS $W$ \+ charm
data on the strange-related distributions, we conduct series of fits to the
HERA I+II $e^{\pm}p$, ATLAS and CMS $W$ \+ charm data together with the
available $W$-lepton charge asymmetry measurement, applying both Eqs. (2) and
(3) using the XFITTER program. After clarifying the reason behind the
contradictory results between ATLAS and CMS data, it is necessary to pin down
what the distribution of strange quark PDFs looks like. To achieve this, we
consider that the combined usage of all the data available here in a fit with
a flexible parametrization form to accommodate the experimental uncertainties
is much preferable to the analysis of merely ATLAS or CMS data using Eqs. (2)
and (3). In this way, the strange quark PDFs could be even more constrained.
Thus, we do one more round of fit with the following parametrization form at
$Q_{0}^{2}=1.9~{}{\mathrm{GeV^{2}}}$ including all the data at once

$\displaystyle xg(x)$
$\displaystyle=A_{g}x^{B_{g}}(1-x)^{C_{g}}-A^{\prime}_{g}x^{B^{\prime}_{g}}(1-x)^{C^{\prime}_{g}},$
(4a) $\displaystyle xu_{v}(x)$
$\displaystyle=A_{u_{v}}x^{B_{u_{v}}}(1-x)^{C_{u_{v}}}(1+E_{u_{v}}x^{2}),$
(4b) $\displaystyle xd_{v}(x)$
$\displaystyle=A_{d_{v}}x^{B_{d_{v}}}(1-x)^{C_{d_{v}}},$ (4c) $\displaystyle
x\bar{u}(x)$
$\displaystyle=A_{\bar{u}}x^{B_{\bar{u}}}(1-x)^{C_{\bar{u}}}(1+D_{\bar{u}}x),$
(4d) $\displaystyle x\bar{d}(x)$
$\displaystyle=A_{\bar{d}}x^{B_{\bar{d}}}(1-x)^{C_{\bar{d}}}(1+D_{\bar{d}}x),$
(4e) $\displaystyle x\bar{s}(x)$
$\displaystyle=xs(x)=A_{\bar{s}}x^{B_{\bar{s}}}(1-x)^{C_{\bar{s}}}(1+D_{\bar{s}}x).$
(4f)

In this parametrization form, we add an extra parameter $D$ to the sea quarks
and put no restrictions between the individual PDF parameters, so that it
might adapt the tensions between individual partons and measurements well.
Setting $s=\bar{s}$ is the usual way used by many PDFs working groups because
of the limited data constraints on the strange quark distribution. By taking
difference of $\sigma(W^{+}\bar{c})$ and $\sigma(W^{-}c)$, it is possible to
extract information on $x(s-\bar{s})$. Recently, the newly extracted PDFs
MSHT20 [18] and the reference [27] obtain very small, order of $10^{-3}$,
$x(s-\bar{s})$ distribution by analyzing both the dimuon and LHC data. Given
some limitations in data on the strange quark distribution, we use symmetric
strangeness throughout in our analysis. Actually, it was predicted [44] that
the momentum and helicity distributions of the strange and antistrange quarks
are different from each other. Hopefully, this problem will be solved by
analyzing relevant data from the near-future high-luminosity and high-energy
experiments and with more precise theoretical calculations.

The denomination for each round of fit goes as follows:

  * •

Fit1.—Applying the parametrization form of Eq. (2) with HERA I+II $e^{\pm}p$
[14] and ATLAS $W/Z$ data [15].

  * •

Fit2.—Applying the parametrization form of Eq. (3) with HERA I+II $e^{\pm}p$
[14] and CMS $W$ \+ charm and $W$-lepton charge asymmetry data [20, 21, 22,
23].

  * •

Fit1p.— Adding ATLAS $W$ \+ charm data [19] to Fit1.

  * •

PKU-fit.—Applying Eq. (4) with all the data.

To get a data-driven characteristics of the PDFs, we get rid of all the
restrictions set on the fit parameters, $A_{\bar{u}}$ = $A_{\bar{d}}$ and
$B_{\bar{u}}$ = $B_{\bar{d}}$ = $B_{\bar{s}}$ in the above fits. The initial
factorization and renormalization scale is still 1.9 $\mathrm{GeV^{2}}$. The
heavy quark masses are set $m_{c}$ = 1.43 GeV and $m_{b}$ = 4.5 GeV. The
strong coupling constant $\alpha_{s}$ is set to be $0.111$ suggested in the
HERAPDF2.0 [14] through a $\chi^{2}$ scan on $\alpha_{s}$. We did conduct a
few rounds of fits with $\alpha_{s}$ values equal to 0.110, 0.111, 0.113,
0.115, 0.117, 0.118 using Eq. (2), imposing specific cuts on $Q^{2}\gtrsim$
3.5, 7.5, 10 $\mathrm{GeV^{2}}$, and found that fit results are almost
independent on the strong coupling constant. To use as much as data, we set
minimum $Q_{\mathrm{min}}^{2}$ = 3.5 $\mathrm{GeV^{2}}$ as a universal cut on
the squared 4- momentum transferred. The Quality of the fit is judged by
$\chi^{2}$/DOF (degrees of freedom (DOF), defined as the difference of the
number of data points used between the number of free fit parameters), where
$\chi^{2}$ is constructed by including the measurement, corresponding
theoretical calculation, and the measurement uncertainties (statistical,
correlated and uncorrelated systematic). The detailed form of the $\chi^{2}$
function is constructed as

$\begin{split}\chi^{2}(M,\lambda)&=\sum_{i=1}^{N_{\mathrm{data}}}\frac{1}{\Delta_{i}^{2}}\left(M_{i}+\sum_{\alpha=1}^{N_{\mathrm{sys.corr}}}\Gamma_{\alpha
i}T_{i}\lambda_{\alpha}-T_{i}\right)^{2}+\sum_{\alpha=1}^{N_{\mathrm{sys.corr}}}\lambda_{\alpha}+\\\
&\sum_{i=1}^{N_{\mathrm{data}}}\mathrm{ln}\frac{\Delta_{i}^{2}}{M_{i}^{2}(\delta_{i,\mathrm{stat}}^{2}+\delta_{i,\mathrm{uncor}}^{2})},\end{split}$
(5)

where $M_{i}$ and $T_{i}$ stand for the measured and theoretical values of the
cross sections. $\Gamma_{\alpha i}$, together with the nuisance parameter
$\lambda_{\alpha}$, quantifies the contribution of each correlated systematic
error source $\alpha$. $\delta_{i,\mathrm{stat}}$ and
$\delta_{i,\mathrm{uncor}}$ are the corresponding relative statistical and
uncorrelated systematical errors which are proportional to measured value. The
index $i$ runs over all $N_{\mathrm{data}}$ data points while index $\alpha$
runs all of the correlated error sources. Because of the multiplicative nature
of both the statistical and uncorrelated systematical errors, $\Delta_{i}^{2}$
is defined as

$\Delta_{i}^{2}=\delta_{i,\mathrm{stat}}^{2}T_{i}M_{i}+\delta_{i,\mathrm{uncor}}^{2}M_{i}^{2},$
(6)

which is corrected for statistical fluctuations in data by scaling with
predicted values according to the recipe in Ref. [45]. This form of
$\Delta_{i}^{2}$ leads to a logarithmic term arising from likelihood
transition of $\chi^{2}$. The minimization of $\chi^{2}$ function against the
fit parameters and the nuisance parameter is done using the package MINUIT
[46], which is interfaced with XFITTER. The PDFs uncertainties arising from
the uncertainties in the measurements are estimated using Hessian method [47]
adopting the same tolerance $\Delta\chi^{2}$ = 1 as the HERA, ATLAS and CMS
groups, corresponding to 68% confidence level (CL).

##  IV Results

Table 1 lists the corresponding $\chi^{2}$/DOF values for four fits showing
the datasets are reasonably agreeable.

Table 1: The corresponding partial, total $\chi^{2}$/DOF function for the four fits Dataset | Fit1 | Fit2 | Fit1p | PKU-fit
---|---|---|---|---
HERA I+II CC $\mathrm{e^{+}p}$ | 49 / 39 | 49 / 39 | 47 / 39 | 44 / 39
HERA I+II CC $\mathrm{e^{-}p}$ | 64 / 42 | 55 / 42 | 66 / 42 | 68 /42
HERA I+II NC $\mathrm{e^{+}p}$ | 226 / 159 | 222 / 159 | 226 / 159 | 226 / 159
HERA I+II NC $\mathrm{e^{+}p,E_{p}}=820$ GeV | 66 / 70 | 70 / 70 | 66 / 70 | 66 / 70
HERA I+II NC $\mathrm{e^{+}p,E_{p}}=920$ GeV | 451 / 377 | 449 / 377 | 450 / 377 | 441 / 377
HERA I+II NC $\mathrm{e^{+}p,E_{p}}=460$ GeV | 220 / 204 | 219 / 204 | 221 / 204 | 220 / 204
HERA I+II NC $\mathrm{e^{+}p,E_{p}}=575$ GeV | 222 / 254 | 224 / 254 | 222 / 254 | 220 / 254
CMS W-$\mu$ charge asymmetry $\sqrt{s}=7$ TeV |  | 14 / 11 |  | 13 / 11
CMS W-$\mu$ charge asymmetry $\sqrt{s}=8$ TeV |  | 3.2 / 11 |  | 5.5 / 11
CMS W + charm $\sqrt{s}=7$ TeV |  | 2.6 / 5 |  | 1.5 /5
CMS W + charm $\sqrt{s}=13$ TeV |  | 2.5 / 5 |  | 1.2 / 5
ATLAS $Z/\gamma^{*}$, $46<m_{\mathrm{z}}<66$ GeV | 23 / 6 |  | 23 / 6 | 24 / 6
ATLAS central $Z/\gamma^{*}$, $66<m_{\mathrm{z}}<116$ GeV | 14 / 12 |  | 15 /12 | 14 / 12
ATLAS forward $Z/\gamma^{*}$, $66<m_{\mathrm{z}}<116$ GeV | 5.6 / 9 |  | 6 / 9 | 4.7 / 9
ATLAS central $Z/\gamma^{*}$, $116<m_{\mathrm{z}}150$ GeV | 6.9 / 6 |  | 6.7 / 6 | 6.8 / 6
ATLAS forward $Z/\gamma^{*}$, $116<m_{\mathrm{z}}150$ GeV | 4.5 / 6 |  | 4.4 / 6 | 4.1 / 6
ATLAS $W^{-}\to l^{-}\bar{\nu}$ | 8.7 / 11 |  | 10 / 11 | 9.3 / 11
ATLAS $W^{+}\to l^{+}\nu$ | 13 / 11 |  | 13 / 11 | 14 / 11
ATLAS $W^{+}$ \+ $\bar{c}$-jet $\sqrt{s}=7$ TeV |  |  | 7.5 / 11 | 7.1 / 11
ATLAS $W^{-}$ \+ $c$-jet $\sqrt{s}=7$ TeV |  |  | 1.6 / 11 | 1.8 / 11
Correlated $\chi^{2}$ | 117 | 90 | 120 | 118
Log penalty $\chi^{2}$ | -14.18 | -9.35 | -14.13 | -6.14
Total $\chi^{2}$ / DOF | 1476 / 1189 | 1391 / 1161 | 1491 / 1211 | 1503 / 1239

Figure 1: The strange PDFs and strange-to-nonstrange ratio as functions of $x$
evaluated at the factorization scale $Q^{2}=1.9~{}\mathrm{~{}GeV^{2}}$ and
$Q^{2}$ = $m_{W}^{2}$ from Fit1 and Fit2. The upper panel shows the central
distribution with error bands corresponding to 68% CL while the bottom panel
shows the ratio to Fit1.

Figure 2: The strange PDFs and strange-to-nonstrange ratio as functions of $x$
evaluated at the factorization scale $Q^{2}=1.9~{}\mathrm{~{}GeV^{2}}$ and
$Q^{2}$ = $m_{W}^{2}$ from Fit1p (blue hatched band) compared to Fit1 and
Fit2. The upper panel shows the central distribution with error bands at 68%
CL while the bottom panel shows the ratio to Fit1.

Figure 3: The $xs(x,Q^{2})$ and $R_{s}$ distributions evaluated at $Q^{2}$ = 4
$\mathrm{GeV}^{2}$ as functions of $x$ for the PKU-fit in comparison with the
previous results. The bands represent 68% CL of experimental uncertainties.

Figure 4: The measured differential cross section of ATLAS $W$ \+ charm
production as a function of lepton $|\eta|$ compared to predictions computed
with the resulted PDF sets: $W^{-}$+c-jet (left panel) and
$W^{+}$+$\bar{c}$-jet (right panel). The measurements are shown with black
error bars for systematic + statistical uncertainty and with filled bands for
total uncertainty (including uncorrelated certainty). The predictions are
obtained with PDF uncertainties at 68%CL.

Figure 5: The measured differential cross section of CMS $W$ \+ charm
production as a function of lepton $|\eta|$ compared to predictions computed
with the resulted PDF sets:$\sqrt{s}=$ 7 TeV (left panel) and $\sqrt{s}=$13
TeV (right panel). The measurements are shown with black error bars for
systematic + statistical uncertainty and with filled bands for total
uncertainty (including correlated uncertainty). The predictions are obtained
with PDF uncertainties at 68%CL.

In the first place, we show the resulted strange quark PDFs and the ratio
$R_{s}$, with ratio to Fit1, extracted from ATLAS $W/Z$ and CMS $W$ \+ charm
data in Fig. 1 applying Eq. (2) (Fit1) and Eq. (3) (Fit2) at the factorization
scale $Q^{2}$= 1.9 $\mathrm{GeV^{2}}$ and the scale at which the $W$ boson is
produced. The $xs(x,Q^{2})$ distributions display overall agreement within the
uncertainties at both scales. However, it is clear to see that they differ
distinctly from each other toward large $x$. Similarly, the $R_{s}$
distributions on the bottom pictures of Fig. 1 follow the same trend that the
strange quark PDF has, thus indicating that the ATLAS $W/Z$ and CMS $W$ \+
charm data have some tensions at large momentum fraction. The central
distribution of strangeness from the fit of ATLAS $W/Z$ data, Fit1, whether in
terms of parton level or the ratio $R_{s}$, is slightly higher than the one
obtained in Fit2. The ratio $R_{s}$ indeed has a significant $x$-dependent
property that coincides with the results obtained by CMS analysis of $W$ \+
charm data [23]. The incompatibility between strange quarks is around
3$\sigma$.

So, the reason that gives rise to the different interpretations of strangeness
in terms of the ratio $R_{s}$ between the original ATLAS and CMS analyses lies
in the different data structure or physics processes, and, possibly, the
somewhat different parametrization forms with flexible parameters of PDFs also
play some role; if the sea quarks are related with each other by setting
$A_{\bar{u}}$ = $A_{\bar{d}}$ and $B_{\bar{u}}$ = $B_{\bar{d}}$ =
$B_{\bar{s}}$, they shape out the ATLAS-type $R_{s}$ as the ATLAS
Collaboration obtained in [15, 17]. On the other hand, they produce the
similar CMS-type of $R_{s}$ if the parameters are set free. As for the
disparity on this ratio arising toward large momentum fraction in Fit1 and
Fit2, this is mainly because Fit1 does not include the $W$ \+ charm production
data that have more significance than $W/Z$ process when it comes to the
strange quark PDFs.

In the second place, the resulted strange quark PDFs and the ratio $R_{s}$ of
Fit1p are depicted in Fig 2 in comparison with Fit1 and Fit2. After the
inclusion of ATLAS $W$ \+ charm data, $xs(x,Q^{2})$ distribution of Fit1p,
depicted with blue hatched band, demonstrates a decent agreement with Fit1 and
Fit2 at $x\lesssim 10^{-2}$, while being compatible with $xs(x,Q^{2})$ of Fit2
toward large $x$. The $R_{s}$ of Fit1p, shown on the bottom, also manifests
compatibility with that of Fit2 with a sharp deviation from that of Fit1 at
$x\gtrsim 10^{-2}$. This compatibility is highly because that the $W$ \+ charm
data exert more constraints on the strange quark PDFs in Fit1.

Based on these results, we can safely state that the observed disparity in the
interpretations of strangeness by the ATLAS and CMS groups, can be ascribed
mainly to the different usage of data from the different physical processes.
The outcome would be consistent if this factor were handled properly with more
flexible parametrization forms.

The combined usage of the inclusive $W/Z$ measurements with the $W$ \+ charm
data and $W$-lepton charge asymmetry provides a better sensitivity on the
strangeness than used separately. Having clarified the discrepancy and its
origin between ATLAS and CMS interpretations of the strange quark PDFs, a
fairly exact distribution of the strange quark PDFs obtained through the PKU-
fit is shown in Fig. 3 in comparison with Fit1, Fit2, and Fit1p at the initial
and $W$ mass scale. It is clear that the strange quark PDFs of PKU-fit are
still in accordance with that of Fit2 and Fit1p throughout entire $x$ while
agreeing with Fit1 only at $x\lesssim 10^{-2}$ and deviating from it toward
large $x$. As Fit1p and Fit2 display, PKU-fit also manifests a smoothly
$x$-dependent behavior of $R_{s}$ at $Q^{2}=1.9~{}{\mathrm{GeV^{2}}}$, which
is declining faster than other fits toward large $x$. The ratio $R_{s}$
distribution at $Q^{2}$ = $m_{W}^{2}$ tends to unity at $x\lesssim 0.02$,
suggesting $(s+\bar{s})=(\bar{u}+\bar{d})$. Generally, we can confirm that the
strange quark PDFs, from both the ATLAS $W/Z$ and CMS $W$ \+ charm data, are
fairly consistent within uncertainties in smaller momentum fraction regions
and manifest a tension toward large momentum fraction. This tension between
fades away when the same physical processes are used with flexible
parametrization forms. The strange quark density of the proton is only
suppressed at $x>0.023$ and enhanced toward smaller $x$.

As confirmed, the inclusion of $W$ \+ charm data together with $W/Z$ from both
the ATLAS and CMS experiments produce well-constrained strangeness. We
calculate the cross section of $W$ \+ charm production of both the ATLAS and
CMS experiments with the extracted PDF sets to see how much the strange quark
PDFs from PKU-fit improve the theoretical results. This is pictured in Figs. 4
and 5 compared with the measurements. The predictions of both the ATLAS and
CMS $W$ \+ charm measurements agree with the measurements quite well. However,
it is worth noting that the calculation obtained from Fit1 is apparently
bigger than the measurements at $|\eta|\approx 2$. The deviations at the large
pseudorapidity have a relation with strange quark density at the corresponding
momentum fraction $x$. This means that the deviation appears when there is an
increase in the strange quark density. Understandably, it is the large-$x$
distribution of the strange quark that leads to the deviation. This tells us
that the uncertainty of the strange quark density in Fit1 by fitting only the
ATLAS $W/Z$ data is large at large momentum fraction.

##  V Conclusions

The shapes of the strange quark PDFs of the proton extracted from LHC data
were interpreted in different ways: the ATLAS statement of fully enhanced
strangeness and the CMS statement of suppressed one which supports the same
idea obtained from fixed target neutrino-nuclear collision experiments. In
this paper, through four rounds of analyses of the HERA I+II $e^{\pm}p$, LHC
$W/Z$, $W$ \+ charm production, and $W$-lepton charge asymmetry data, i.e.,
Fit1, Fit2, Fit1p and PKU-fit, we have investigated the exact distribution of
the strange quark and its ratio to nonstrange sea quarks implementing the
ATLAS and CMS analyses set up with free sea quark parameters using XFITTER
program, and resolved tensions between ATLAS and CMS data, and finally
extracted fairly well-shaped strange quark densities. At first, both the ATLAS
$W/Z$ and CMS $W$ \+ charm data are proved, through Fit1 and Fit2 fits, to be
partially consistent with the relaxed sea quark parameters in terms of strange
quark and its ratio to nonstrange sea quark distributions at $x\lesssim
10^{-2}$. The inconsistency toward large $x$ between strange quark density and
the ratio $R_{s}$ from the ATLAS and CMS Collaborations suggests that there is
tension between ATLAS 7$~{}\mathrm{TeV}$ $W/Z$ and CMS 7 and
13$~{}\mathrm{TeV}$ $W$ \+ charm data. It is shown in Fit1p that this tension
can be eliminated by fitting both the $W/Z$ and $W$ \+ charm data
simultaneously. The shapes of the strange quark PDFs and the ratio obtained in
our analysis are purely data driven. There is no handmade constraint on the
fit parameters. The $W$ \+ charm cross section is calculated with the obtained
PDFs and deviations from measurement at the large lepton $|\eta|$ are
observed. The PDFs obtained from only the HERA I+II $e^{\pm}p$ and ATLAS $W/Z$
measurement produces larger cross section values than the $W$ \+ charm data at
the large lepton $|\eta|$, suggesting a large uncertainty of the proton
strangeness from fitting only the ATLAS $7~{}\mathrm{TeV}$ $W/Z$ data at large
$x$. When the HERA I+II $e^{\pm}p$, ATLAS $W/Z$, CMS $W$ \+ charm and
$W$-lepton data are used in combination, further strengthened by ATLAS $W$ \+
charm measurement, with a more flexible parametrization form being applied,
analyzed in PKU-fit, the resulted distribution of the strange quark density
gets improved and the ratio $R_{s}$ shows more $x$-dependent property. The
strangeness is gradually suppressed at larger $x$ but it is highly enhanced
toward smaller $x$ with its magnitude becoming comparable to that of the
light-flavor up and down sea quarks.

##  VI Acknowledgment

We thank the XFITTER team for providing useful help and discussions on the
usage of their awesome program. We are also so grateful to Katerina Lipka from
DESY for providing the corresponding grid file of CMS $W$ \+ charm production
at $\sqrt{s}=13~{}\mathrm{TeV}$ for the corresponding theoretical calculation.
This work is supported by the National Natural Science Foundation of China
(Grant No.12075003).

## References

  * [1] E. D. Bloom, D. H. Coward, H. C. DeStaebler, J. Drees, G. Miller, L. W. Mo, R. E. Taylor, M. Breidenbach, J. I. Friedman, G. C. Hartmann and H. W. Kendall, Phys. Rev. Lett. 23, 930 (1969).
  * [2] M. Breidenbach, J. I. Friedman, H. W. Kendall, E. D. Bloom, D. H. Coward, H. C. DeStaebler, J. Drees, L. W. Mo and R. E. Taylor, Phys. Rev. Lett. 23, 935 (1969).
  * [3] L. Del Debbio et al. (NNPDF), J. High Energy Phys. 03, 039 (2007).
  * [4] S. Dulat, T. J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump and C. P. Yuan, Phys. Rev. D 93, 033006 (2016).
  * [5] L. A. Harland-Lang, A. D. Martin, P. Motylinski and R. S. Thorne, Eur. Phys. J. C 75,204 (2015).
  * [6] S. Alekhin, J. Blumlein, S. Klein and S. Moch, Phys. Rev. D 81, 014032 (2010)
  * [7] A. M. Cooper-Sarkar (ZEUS and H1 Collaboration), PoS EPS-HEP2011, 320 (2011).
  * [8] U. Baur, F. Halzen, S. Keller, M. L. Mangano and K. Riesselmann, Phys. Lett. B 318, 544 (1993).
  * [9] E. Farhi and R. L. Jaffe, Phys. Rev. D 30, 2379 (1984).
  * [10] M. Goncharov et al. (NuTeV Collaboration), Phys. Rev. D 64, 112006 (2001).
  * [11] D. Mason et al. (NuTeV Collaboration), Phys. Rev. Lett. 99, 192001 (2007).
  * [12] O. Samoylov et al. (NOMAD Collaboration), Nucl. Phys. B876, 339 (2013).
  * [13] G. Aad et al. (ATLAS Collaboration), Phys. Rev. Lett. 109, 012001 (2012).
  * [14] H. Abramowicz et al. (H1 and ZEUS Collaborations), Eur. Phys. J. C 75, 580 (2015).
  * [15] M. Aaboud et al. (ATLAS Collaboration), Eur. Phys. J. C 77, 367 (2017).
  * [16] A. M. Cooper-Sarkar and K. Wichmann, Phys. Rev. D 98, 014027 (2018).
  * [17] G. Aad et al. (ATLAS), J. High Energy Phys. 07, 223 (2021).
  * [18] S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin and R. S. Thorne, Eur. Phys. J. C 81, 341 (2021).
  * [19] G. Aad et al. (ATLAS Collaboration), J. High Energy Phys. 05, 068 (2014).
  * [20] 20 S. Chatrchyan et al. (CMS Collaboration), Phys. Rev. D 90, no.3, 032004 (2014).
  * [21] V. Khachatryan et al. (CMS Collaboration), Eur. Phys. J. C 76, 469 (2016).
  * [22] S. Chatrchyan et al. (CMS Collaboration), J. High Energy Phys. 02, 013 (2014).
  * [23] A. M. Sirunyan et al. (CMS Collaboration), Eur. Phys. J. C 79, 269 (2019).
  * [24] N. Yalkun and S. Dulat, Chin. Phys. C 43, 123101 (2019).
  * [25] C. Schmidt, J. Pumplin, C. P. Yuan and P. Yuan, Phys. Rev. D 98, 094005 (2018).
  * [26] T. J. Hou, J. Gao, T. J. Hobbs, K. Xie, S. Dulat, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin and C. Schmidt, et al. Phys. Rev. D 103, 014013 (2021).
  * [27] F. Faura, S. Iranipour, E. R. Nocera, J. Rojo and M. Ubiali, Eur. Phys. J. C 80, 1168 (2020).
  * [28] S. Alekhin et al., HERAFitter Eur. Phys. J. C 75, 304 (2015).
  * [29] 29 C. Gong and B.-Q. Ma, Phys. Rev. C 97, 065207 (2018).
  * [30] S. J. Brodsky, I. Schmidt and J. J. Yang, Phys. Rev. D 70, 116003 (2004).
  * [31] V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. 15, 438 (1972).
  * [32] V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. 15, 675 (1972).
  * [33] Y. L. Dokshitzer, Sov. Phys. JETP 46, 641 (1977).
  * [34] G. Altarelli and G. Parisi, Nucl. Phys. B 126, 298 (1977).
  * [35] S. Fanchiotti, B. A. Kniehl and A. Sirlin, Phys. Rev. D 48, 307 (1993).
  * [36] M. Botje, Comput. Phys. Commun. 182, 490 (2011).
  * [37] R. S. Thorne and R. G. Roberts, Phys. Rev. D 57, 6871 (1998).
  * [38] R. S. Thorne, Phys. Rev. D 73, 054019 (2006).
  * [39] R. S. Thorne, Phys. Rev. D 86, 074017 (2012).
  * [40] J. M. Campbell and R. K. Ellis, Phys. Rev. D 60, 113006 (1999).
  * [41] J. M. Campbell, R. K. Ellis and C. Williams, J. High Energy Phys. 07, 018 (2011).
  * [42] T. Carli, D. Clements, A. Cooper-Sarkar, C. Gwenlan, G. P. Salam, F. Siegert, P. Starovoitov and M. Sutton, Eur. Phys. J. C 66, 503-524 (2010).
  * [43] 43 S. Alekhin, J. Blümlein and S. Moch, Eur. Phys. J. C 78 477 (2018)
  * [44] S. J. Brodsky and B.-Q. Ma, Phys. Lett. B 381, 317-324 (1996)
  * [45] F. D. Aaron et al. (H1 and ZEUS Collaborations), J. High Energy Phys. 01 (2010), 109
  * [46] F. James and M. Roos, Comput. Phys. Commun. 10, 343-367 (1975)
  * [47] J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai and W. K. Tung, Phys. Rev. D 65, 014013 (2001).