Adrian Alan Pol
email: adrian.pol@cern.ch

Gianluca Cerminara

Cécile Germain

Maurizio Pierini

Agrima Seth

Cecile Germain

Cecile Germain

Detector monitoring with artificial neural networks at the CMS experiment at the CERN Large Hadron Collider

Keywords: High Energy Physics, Large Hadron Collider, Compact Muon Solenoid, Machine Learning, Data Quality Monitoring, Artificial Neural Networks

Introduction

The Compact Muon Solenoid (CMS) experiment is a general purpose particle physics detector operating at the CERN Large Hadron Collider [START_REF] Germain | The large hadron collider, conceptual design[END_REF] (LHC). Data collected with the CMS detector are used in many aspects of modern particle physics, notably the discovery [START_REF] Chatrchyan | Observation of a new boson at a mass of 125 gev with the cms experiment at the lhc[END_REF] and characterization [START_REF] Vardan Khachatryan | Precise determination of the mass of the Higgs boson and tests of compatibility of its couplings with the standard model predictions using proton collisions at 7 and 8 TeV[END_REF] of the Higgs boson.

The CMS detector is described in details in [START_REF] Chatrchyan | The cms experiment at the cern lhc[END_REF], together with a definition of the used coordinate system and the relevant kinematic variables. In CMS, muons are measured with detection planes instrumented with four detector technologies: drift tubes (DTs), cathode strip chambers, resistive plate chambers, and gas electron multipliers. A detailed description of the CMS muon detectors can be found in [START_REF] Albert M Sirunyan | Performance of the CMS muon detector and muon reconstruction with proton-proton collisions at √ s = 13 TeV[END_REF].

Within the CMS Collaboration, physics analysis are performed on good data, selected by imposing stringent quality criteria. During data taking, a subset of the collected statistics is processed in real time, to create a set of histograms filled with a certain critical quantities. Statistical tests are performed to compare these histograms to a set of predefined reference, representing the typical detector response during normal operation conditions. Using the histogram comparison and the outcome of the tests, expert shifters acknowledge the alarms and may decide to intervene (up to stopping the data taking), depending on the evaluation of the problem severity. The knowledge of the LHC running conditions and of the history of possible issues identified in the past, are key ingredients in this decision process. Details on the infrastructure used for this Data Quality Monitoring (DQM) are given in [START_REF] De | The data quality monitoring challenge at cms: experience from first collisions and future plans[END_REF]. The two main domains of the monitoring chain are:

online monitoring, which provides live feedback on the quality of the data while they are being acquired, allowing the operator crew to react to unforeseen issues identified by the monitoring application;

offline monitoring, designed to certify the quality of the data collected and stored on disk using centralized processing (referred to as the event reconstruction, which converts detector hits into a list of detected particles, each associated with an energy and direction).

These two validation steps differ in three main aspects:

the latency of the evaluation process; online monitoring is requested to identify anomalies in quasi real time to allow the operators to intervene promptly while the offline procedure has a typical timescale of several days, the fraction of the data which they have access to;

online processing runs at a rate of 100 Hz, corresponding to approximately 0.1% of the data written to disk for analysis, while the offline processing takes as input the full set of events accepted by the trigger system (∼ 1 kHz of data), the granularity of the monitored detector components; while offline monitoring requires identifying only overall status of the detector components, online should determine faulty subdetector elements.

Despite their specific characteristics, these two steps rely on the same anomaly detection strategy: the scrutiny of a long list of predefined histograms, selected to detect a set of known failure modes. These histograms are monitored by detector experts, who compare each distribution to a corresponding reference, derived from good data in line with predetermined validation guidelines. This two-layer monitoring protocol was adopted by the CMS Collaboration for LHC Run I (2010-2012) and in Run II (2015)(2016)(2017)(2018). The ever increasing detector complexity, monitoring data volumes and the necessity to cope with different LHC running scenarios call for an increasing level of automation of the applications in the future. Already, the amount of histograms to monitor is challenging for a single shifter, while the number of histograms to monitor increases every time a new failure mode is identified and consequently added to the list of known potential problems. Furthermore, the human intervention and currently implemented tests require collecting a substantial amount of data, implying a detection delay. Last but not least, the cost in terms of human resources is substantial i.e. the 24/7 DQM shifter and the expert personnel responsible for updating the good data references and related instructions. We believe that introducing machine learning into the CMS DQM process will help with those challenges.

This work focuses on the online monitoring. We concentrate on the application of deep learning techniques, and specifically image-like processing [START_REF] Lecun | Deep learning[END_REF] for the au-tomation of detector level monitoring. While the main focus of this work is on improving detection specificity and sensitivity, the proposed approach could come with practical advantages in operation being based on less astringent assumption on the nature of the anomalies.

As a concrete example we use real data recorded by the CMS DT chambers of the muon spectrometer during the data-taking campaign of the LHC Run II. The main aspects of this work are:

we exploit the geographical information of the detector assessing the (mis)behavior with high-granularity and then combining the results to probe different detector components; we detect different types of anomalies affecting the detector at different scales (ranging from a few channels to collective behaviors of big portion of the DT system) by combining different algorithms; we show that image-like processing achieves considerably better performance with respect to the current threshold-based DQM test (later called the production algorithm) and allows to tune the working point in terms of specificity (depending on the deployment strategy).

Although the experimental demonstration of the results presented in this paper is tied to the specificities of the DT subdetector, the procedure that we discuss have a potentially broader application scope. Mainly because the typical issues encountered with other subdetectors are analogous. This possibility is currently under investigation for other detector components of the CMS experiment.

The remainder of the paper is organized as follows. Sections 2 and 3 present in more details the problem that we want to solve and describe the utilized data set. Section 4 reviews the current state of the art in the machine learning domain of failure detection. Sections 5, 6 and 7 present three complementary approaches to the problem. Section 8 describes and discusses the results.

The CMS Drift Tube muon system

An illustration of the internal structure of a DT chamber is shown in Fig. 1. Each chamber, on average 2 × 2.5 m in size, consists of 12 layers of drift tubes. Layers are arranged in three groups of four, each containing a variable number of tubes, up to 96. The middle group measures the coordinate along the direction parallel to the beam and the two outer groups measure the perpendicular coordinate. Each tube corresponds to one readout channel (briefly referred to as channel in the rest of the paper). Particles carrying an electromagnetic charge and traversing a tube release an electronic signal by ionizing the gas in the tube (a hit). By combining the information provided by the channels, one can determine the trajectory of the particle crossing the chamber.

The chamber numbering schema follows that of the iron of the yoke, consisting in five wheels (see Fig. 2) along the z-axis, each one divided into 12 azimuthal sectors (see Fig. 3). The wheels are numbered from -2 to +2, sorted according to global CMS z-axis, with wheel 0 situated in the central region around the proton-proton collision point. The sector numbering is assigned in an anti-clockwise sense when looking at the detector from the positive z-axis, starting from the vertically-oriented sector on the positive-x side in the CMS coordinate system (sector 1). Chambers are arranged in four stations at different radii, named MB1, MB2, MB3 and MB4. The first and the fourth stations are mounted on the inner and outer face of the yoke respectively; the remaining two are located in slots within the iron. Each station consists of 12 chambers (one per sector) except for MB4 (which contains 14 chambers). The total number of chambers is then 5 × (3 × 12 + 14) = 250. 3 The data set, monitoring strategy, and preprocessing 3.1 The occupancy matrix CMS data are organized in acquisition runs (or just runs in CMS jargon), corresponding to a given setup both of the CMS detector and of LHC accelerator. Runs are denoted by integers, increasing with time. Their duration is varying from as little as few seconds to as much as several hours.

Each run is divided into luminosity sections (LSs), a time interval corresponding to a fixed number of protonbeam orbits in the LHC and amounting to approximately 23 seconds, numbered progressively from 1 at the start of each run. Each LS can be identified uniquely by specifying the LS number and the run number. The beam intensity (also referred to as luminosity) varies along each run, resulting in a varying number of protonproton collision data (the events).

For each chamber k in a given run, the current DQM infrastructure [START_REF] Tuura | Cms data quality monitoring web service[END_REF] records an occupancy matrix C k , which contains the total number of particle hits at each channel for a given LS or set of consecutive LSs. The occupancy matrix can be viewed as a varying size two-dimensional array organized along layer (row) and channel (column) indices:

C k = {x k i,j ; 1 ≤ i ≤ l, 0 ≤ j < n i },
where l = 12 is the number of layers and n i is the number of channels in layer i. In general, we label the chambers and their components as C k and x k i,j . For simplicity, we omit the k index when discussing problems related to individual chambers, until Section 6. Figure 4 shows examples of occupancy matrices, represented as two-dimensional occupancy plots. The utilized data set consists of 21000 occupancy matrices for the 250 chambers.

Monitoring strategy

The anomaly detection method currently used in the online monitoring production system targets a specific failure scenario: a region of cells not providing any electronic signal, large enough to affect the track reconstruction in the chamber. This is by far the most frequent issue, usually related to transient problems in the readout electronics. Examples of this kind of failures are shown in Fig. 4 B andC. These kinds of occupancy plots are created accumulating data in time. Once in a while, the plot filling process is reset, to increase sensitivity to problems occurring during the run. The production algorithm evaluates samples per chamber. Although it quantifies the fault severity on the basis of the fraction of affected channels, it does not identify specific faulty layers.

The novel approach proposed in this work goes beyond the functionalities of the current production algorithm. Starting from the identification of layers with under-performing cells, it provides effective identifica-tion of faulty chambers. Moreover, it exploits the geographical information of the layer and chamber position to identify different kind of failures. To this purpose, three complementary approaches are considered:

-Local: data collected in each layer are treated independently from the others. As for the production algorithm, this approach regards chambers which have occupancy of hits with small variance between neighboring channels as expected behavior and targets a well known list of problems with a supervised approach. Chambers which have dead, inefficient or noisy regions, are considered problematic. We explore this approach in Section 5.

-Regional: we extend the local approach to account for intra-chamber problems, to be applied whenever faults are spotted only when the information about all layers within one chamber is present. For this purpose, we simultaneously consider all layers in a chamber, but each chamber is considered independently from the others, Section 6.

-Global: we simultaneously use the information of all the chambers for a given run. The position of the chamber in the CMS detector (uniquely determined by the wheel, station, and sector numbers) impacts expected occupancy distribution of the channel hits. This approach is described in Section 7.

Preprocessing

A common data set preprocessing procedure is used for the three studies (for visual interpretation, see Fig. 5).

-Standardization of the chamber data: the number of channels x in a layer varies not only within the chamber but also depends on the chamber position in the detector. This quantity falls between 47 and 96. We enforce a fixed-input dimensionality by applying a layer-by-layer one dimensional linear interpolation to match the size of the smallest layer s in data set. The smallest layer is chosen to simplify our models later in this study. Starting from the recorded matrix x ij , a standardized matrix xi,j is defined as:

xi,j = f rac(α)(x i, α -x i, α ) + x i, α .
where α is an interpolation point, defined by α = j ni ns . We verified that this method doesn't compromise sensitivity to very small problematic regions despite a small reduction in the amplitude and sharpness of the anomalies.

-Smoothing: according to CMS DT experts, misbehaving channels are problematic only when a spatially contiguous cluster of them is observed. Instead, isolated misbehaving channels are not considered a problem. To take this into account the one dimensional median filter is applied:

xi,j = med(x i,j , x i,j+1 , x i,j+2 ).

-Normalization: the occupancy of the chambers in the input data set depends on the integration time and on the LHC beam configuration and intensity i.e. on the number of LSs spanned when creating the image and corresponding luminosity. The normalization strategy depends on the need for comparing data across chambers or across runs: the precise procedure used in the two approaches is described in Sections 5 and 6, respectively.

Machine learning for DQM anomaly detection

In this section we briefly discuss machine learning anomaly detection techniques in light of both the operational condition and the a priori knowledge of the data. Machine learning presents several advantages over the currently adopted procedure as the decision function can be learned from collected data. In the future, it might be possible to bypass human intervention when the algorithm decision is not controversial and only invoke the shifters opinion for intermediate questionable cases. An example of this approach is discussed in [START_REF] Borisyak | Towards automation of data quality system for cern cms experiment[END_REF] in the context of the CMS offline monitoring. The high data dimensionality precludes simple parametric density estimation of the normal behavior. This leaves an extremely wide range of methods such as one-class Support Vector Machine (µ-SVM) [START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF], Isolation Forest [START_REF] Fei | Isolation forest. In Data Mining[END_REF][START_REF] Fei | Isolation-based anomaly detection[END_REF] and different flavors of deep learning. For a general survey see [START_REF] Charu | Outlier analysis[END_REF].

Anomaly detection techniques usually assume rarity of abnormal events (considered as outliers with respect to the normal generating process) and/or lack of a complete set of representative examples of all possible behaviors. If such representative examples are available, anomaly detection reduces to binary classification (supervised learning), with possibly the help of various resampling methods [START_REF] Charu | Data classification: algorithms and applications[END_REF] or reformulation of the objective function [START_REF] Cowan | Asymptotic formulae for likelihood-based tests of new physics[END_REF] for dealing with class imbalance.

In our case, supervised learning is clearly a valid option as specific anomalous scenarios were extensively studied. The CMS DQM framework keeps copious archives of subdetector-specific quality-related quantities, e.g. the DT occupancy plots. Moreover, the imbalance between good and bad data is not extreme, with a typical rate of anomalies reaching the 10% level. These anomalies are then frequent enough for a sizable set of them to be used for supervised training.

However, there is a good motivation for a semisupervised anomaly detection approach. Beside the deep learning methods that will be discussed at the end of this section, we experiment with the two reference methods, which are variants of one-class classification: µ-SVM and Isolation Forest. µ-SVM estimates the support of the data distribution by a non-linear (kernel) transform of the data space (as in all SVM techniques) and by identifying the hyperplane that maximizes the separation of the training data from the origin. Accordingly, µ-SVM has the important property of being a novelty detection algorithm: once trained, it is not sensitive to the frequency of anomalies. However, the implicit prior of kernel-based classification is that the function to be learned is smooth such that generalization can be achieved by local interpolation between neighboring training examples. As argued at length by Bengio et al. (for instance in [START_REF] Bengio | Scaling learning algorithms towards ai[END_REF] and [START_REF] Bengio | Representation learning: A review and new perspectives[END_REF]), this assumption is questionable for high data dimensionality. An alternative is the Isolation Forest, which copes with the curse of dimensionality by relying only on the principle of isolation of outliers in a random recursive partitioning of the feature space along the axes and tree ensembles. The Isolation Forest algorithm does not rely on any distance or density measure, but assumes that anomalies can be isolated in the native feature space. Besides being highly scalable to large data sets, Isolation Forest offers some possibility of interpretation.

Classical fully unsupervised approaches based on neighborhood (e.g. k nearest neighbors), topological density estimation (Local Outlier Factor and its variants) or clustering (for a detailed presentation, see [START_REF] Goldstein | A comparative evaluation of unsupervised anomaly detection algorithms for multivariate data[END_REF]) are not relevant here. These algorithms have quadratic complexity and poorly perform in high dimensions, because of data sparsity (in high dimensions, all pairs of points become almost equidistant) [START_REF] Zimek | A survey on unsupervised outlier detection in highdimensional numerical data[END_REF]. Moreover, a simple geometric (e.g. Euclidean) distance in the feature space does not define a similarity metric. For instance, the distance between examples A and B in Fig. 4 is dominated by the contribution of well-behaving channels. The similarity function, or equivalently the adequate representation, must be learned from the data.

This representation learning view [START_REF] Bengio | Representation learning: A review and new perspectives[END_REF] points towards deep learning, as it should remain sensitive to the local geometric relationship in the data related to the underlying apparatus. Convolutional networks [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] integrate the basic knowledge of merely the topological structure of the input dimensions and learn the optimal filters that minimize the objective error.

A more ambitious goal is to extract an explanatory representation of the anomalies with latent variables, in a probabilistic framework (e.g. restricted Boltzmann machines, or variational autoencoders [START_REF] Jimenez Rezende | Stochastic backpropagation and approximate inference in deep generative models[END_REF]), where the learned representation is the posterior distribution of the latent variables given an observed input. However, even the inference step with these representations may suffer from high computational cost, and requires some further feature construction.

The alternative is the trade-off between interpretability and simplicity by learning a direct encoding, typically as a neural network based autoencoder, which is a parametric map from inputs to their representation. Although it has been argued that, even for basic neural networks, most of the training is devoted to learning a compressed representation [START_REF] Tishby | Deep learning and the information bottleneck principle[END_REF][START_REF] Shwartz | Opening the black box of deep neural networks via information[END_REF], autoen-coders are particularly suitable for anomaly detection. When trained on the inliers, testing on unseen faulty sample tend to yield sub-optimal representations, indicating that a sample is likely generated by a different process. In order to go beyond simple dimensionality reduction while preventing over-fitting, various flavors of regularization are proposed (the literature being considerable, we give only some entry points):

sparse autoencoders [START_REF] Marc | Efficient learning of sparse representations with an energy-based model[END_REF] penalize the output of the hidden unit activations or the bias; denoising autoencoders [START_REF] Vincent | Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion[END_REF] robustify the mapping by requiring it to be insensitive to small random perturbations; contractive autoencoders [START_REF] Rifai | Contractive auto-encoders: Explicit invariance during feature extraction[END_REF] pursue the same goal, by penalizing analytically the sensitivity of learned features in a data-driven interpretation of the Tangent Propagation algorithm [START_REF] Simard | Transformation invariance in pattern recognition -tangent distance and tangent propagation[END_REF].

In fact, denoising and contractive autoencoders learn density models implicitly, through the estimation of statistics or through a generative procedure [START_REF] Alain | What regularized auto-encoders learn from the data-generating distribution[END_REF].

5 Local approach: detecting faulty behavior within a layer

Motivation

The first experiment concentrates on training a classifier to identify local problems, i.e. considering each layer independently from the others. This approach enforce the expert knowledge of what is currently considered correct or anomalous and probes the detector with higher granularity than the production algorithm.

The goal is to identify regions of channels not registering any hits (called dead channels in detector jargon), or having lower detection efficiency (hence lower hit counts with respect to the neighboring ones in the same layer) or being dominated by electronic noise (called noisy channels). These are by far the most frequent failure modes. The local approach can be considered as an initial benchmark comparing fully supervised, semi-supervised and unsupervised methods, and specific algorithms in each category, before embarking in full-fledged anomaly detection. Moreover, the local approach, if successful, can be further exploited as a preprocessing step for filtering these trivial faults before attempting to detect more elusive ones.

Given the locality restriction of this approach, contextual information is not accessible. As a consequence of this, a model based on this strategy will not be able to spot, for example, a faulty layer in which occupancy is decreased uniformly with respect to neighboring lay-ers. We acknowledge this limitation and address it in Section 6.

Data set and methods

After having applied the standardization procedure (see Section 3.3), a layer is represented as a single row of an occupancy matrix:

X i = (x i,1 , xi,2 , . . . , xi,47 ).
The available data set consists of 21000 chambers corresponding to 228480 individual layers.

Hit counts in a layer are normalized to a [0, 1] range, dividing them by the maximum of the occupancy value in the layer: ẋi,j = xi,j max(X i ) .

The need for normalization comes from the intrinsic variation of the occupancy, which depends on the spatial position of the chamber (as described in more detail in Section 7) and on the integration time of the analyzed image.

In this experiment, we compare the performances of the following:

unsupervised learning with (a) a simple statistical indicator, the variance within the layer, and (b) an image processing technique, namely the maximum value of the vector obtained by applying a variant of an edge detection Sobel filter [START_REF] Sobel | An isotropic 3×3 image gradient operator[END_REF]:

S i = max( -1 0 1 X i ).
semi-supervised learning, with (c) Isolation Forest, and (d) µ-SVM. supervised learning, with (e) a fully connected shallow neural network (SNN), and (f) a convolutional neural network (CNN);

The ground truth is established by field experts on a random subset of the data set, by visually inspecting the input sample before any preprocessing: 5668 layers were labeled as good and 612 as bad. The 9,75% fault rate is a faithful representation of the real problem at hand. With this ratio, both anomaly and outlier detection approach can be considered. Out of this set 1134 good and 123 bad examples are reserved to compose the test set corresponding to 20% of the labeled layers. The remaining examples are used for training and validation for the semi-supervised and supervised methods.

The Isolation Forest and µ-SVM models are crossvalidated using five stratified data set folds to search for their corresponding optimal hyper-parameters. Subsequently, the Isolation Forest is retrained using those hyper-parameters (100 base estimators in the ensemble) on the full unlabeled data set, while µ-SVM (RBF kernel, ν of 0.4, γ of 0.1) is retrained using only negative class examples. The architecture of the CNN model with one dimensional convolution layers used for this problem is shown in Fig. 6. Rectified linear units are chosen as activation functions for inner-layer nodes, while the softmax function is used for the output nodes. The model is trained using the Adam [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF] optimizer and early stopping mechanism monitoring validation set (set to 20% of data set) with patience set to 32 epochs. The model is implemented in Keras [START_REF] Chollet | [END_REF], using TensorFlow [START_REF] Abadi | Tensorflow: Large-scale machine learning on heterogeneous distributed systems[END_REF] as a backend.

The SNN model consists of one hidden fully-connected layer with 16 units (chosen to approximately match number of parameters in the CNN). As for CNN, it uses rectified linear unit as activation function of the hidden nodes and the softmax function is used for the output nodes. This model is primarily introduced to obtain a term of comparison for the CNN.

Unlike what was done for the other models, we do not apply the smoothing preprocessing step described in Section 3.3 for CNN nor SNN models, in order to allow them to learn their filters. Additionally, we weight our negative S 0 and positive S 1 samples to account for class imbalance. The weight λ ψ for a sample in class ψ ∈ {0, 1} is defined by

λ ψ = |S| 2 • |S ψ |
where S = S 0 ∪S 1 . We discuss the results in Section 8.1.

6 Regional approach: detecting unusual behavior within a chamber

Motivation

In normal conditions, healthy chambers show similar occupancy levels in neighboring layers. The four central chambers have a different behavior due to their different orientation (see Section 2). The regional approach exploits the relative occupancy patterns of the layers within a chamber. For example, it aims at detecting failure modes where the occupancy of hits decreases uniformly in a specific layer or set of layers. Typical examples of these kind of failures are problems related with the high-voltage bias of the drift cells. The voltage distribution system is organized by layers and a lower value w.r.t to the nominal operation point would result in lower detector efficiency and, as a consequence, lower absolute occupancy in the affected region. Figure 7 shows an example of such an occurrence, where layer 9 is misbehaving. The production algorithm and the local models of Section 5 are not conceived to detect this type of anomalies.

Data set and methods

In the early stages of this work we observed that a model capable of detecting regional anomalies cannot be successfully trained if the local faults are not filtered beforehand. Moreover, the available ∼ 500 labeled images do not provide a sufficiently large training set. Thus we start with a much larger data set (all the unlabeled samples). We solve the labeling problem using the score of the convolutional model presented in Section 5 as an approximation of the ground truth. For this, we choose a working point with 99% true positive rate (to guarantee a large data set size) and 5% false positive rate. Approximately 90% of the collected data are labeled as good by online and offline monitoring experts. We then estimate the residual bad example contamination to be ∼ 0.5%. We believe that the residual contamination of problematic chambers is reduced to a tolerable level. All chambers with any layer identified as faulty are discarded. Chambers located in MB4 are discarded as well, because of the lack of a middle group of four layers, see Section 2. The above changes effectively narrowed the training data set to 8441 matrices. The smoothing and standardization procedures are applied to all the layers Ck within each chamber obtaining matrices of shape 12 × 46. The occupancy of hits within This normalized values to the [0, 1] range while retaining the information about the relative occupancy between the layers.

In order to evaluate the model, we use the only labeled set for the class of anomalies that we want to tackle: a subset of the data (runs 302634 and 304737-304740), during which layer 9 of some chambers was operating at a voltage lower than the nominal one (see Fig. 7). In particular, the voltage was set to 3450 V in runs 304737-304740 and to 3200 V in run 302634 while the standard operation point is 3550 or 3600 V depending on the chamber. These settings result in an absolute difference in hit counting, more pronounced for the lower voltage settings, because the physics of gas ionization by radiation. The chambers where all layers operate at nominal conditions are considered as good in the test.

In this experiment, the following semi-supervised approaches are considered:

simple bottleneck autoencoder with the representation layer equal to 20 units; convolutional autoencoder; denoising autoencoder in which we add additional artificial noise to training samples; autoencoder with kernel L1 (10 -5 ) sparsity regularization in the hidden layers. Similarly to the local approach we train the autoencoders using the Adam optimizer. Early stopping mechanism with the patience set to 32 epochs is adopted to monitor validation set (20% of the total data set). All models are implemented using the Keras library with TensorFlow as a backend. The architecture of the model is shown in Fig. 8. A and B for, respectively, the convolutional autoencoder and the other three models (for which a common architecture is adopted). The bottleneck architecture is kept for both denoising and sparse autoencoders in order to limit the amount of parameters to train. The parametric rectified linear unit is used as the activation function on the hidden layers, while the output layer uses the sigmoid function. All models are instructed to minimize the mean squared error (MSE) between original, ẋ, and reconstructed, ẍ, samples:

k = 1 ij i,j ( ẋk i,j -ẍk i,j ) 2 .
We discuss the results in Section 8.2.

7 Global approach: detecting unusual behavior using global information

Motivation

In the third approach, we aim at detecting anomalies looking at the global ensemble of muon chambers, exploiting the dependency of the occupancy of each of them on their position in the detector. We categorize the chambers according to their position in the spectrometer and its impact on the occupancy pattern, exploiting the field knowledge to predetermine the classes of chambers.

The expected occupancy pattern is mainly driven by the proximity to the beam-collision point, at the center of the detector. As a consequence, chambers in different stations (see Section 2) will manifest a different behavior. The rotational symmetry of the detector geometry and of the collision events around the beam axis is taken into account grouping chambers within the same station, independently on the sector they belong to (see Section 2). Similarly chambers belonging to the same station but in opposite wheels are considered alike. Additionally, the behavior of the chambers is expected to be the same across different runs, modulo the overall decrease of occupancy due to the decrease of beam intensity across the fill. This leaves us with a categorization based on the chamber numbering schema, where the station number and the absolute value of the wheel number are the only relevant parameters.

The problem is clearly contextual, in the sense that important explanatory attributes are not part of the basic data features. Conditional anomaly detection [START_REF] Song | Conditional anomaly detection[END_REF] has been proposed to deal with such a situation when the relevance of external attributes is unknown. For instance, if a set of environmental or technical attributes are monitored that can impact the behavior of the detector components. In our case, the spatial position of the chambers are our only external attribute, and their impact is assured by common understanding of the underlying physics processes. Thus, we are back to a point anomaly problem.

Methods

In this approach we use a bottleneck autoencoder similar to that introduced in Section 6 (see Fig. 8), except that the size of the bottleneck layer is reduced to three units for visualization purposes. We also follow the same preprocessing, training and validation procedure. The goal of the study is to exploit the categorization of the chambers based on their geographical location to interpret the compressed representation of the network.

Global faults are not tracked before by DT experts. In absence of a global label, we only considered an unsupervised method for this experiment. We discuss the results in Section 8.3.

Results and Discussion

Local approach

The performance of the various models on a held out test data set can be seen in Fig. 9, where we show the different receiver operating characteristic (ROC) curve. Compared to statistical, image processing or other machine learning based solutions, supervised deep learning clearly outperforms the rest. Thanks to the limited number of parameters of the model, the training converges to a satisfactory result (Fig. 10), despite the number of training samples being small.

Although the Area Under Curve (AUC) of the fullyconnected shallow neural network is comparable to the one of CNN, the latter is a better solution when requiring maximum specificity (true negative rate (TNR), aims at avoiding false positives) and sensitivity (true positive rate (TPR), aims at avoiding false negatives). The relatively good performance of the basic and unsupervised variance method, compared to the poor results of the filter, and the near optimal performance of the SNN, show that the features to learn are not simple contrasts, although the superior performance of the CNN demonstrate that the initial edge detection layer is useful. The limited performance of Isolation Forest is likely to come from the violation of its fundamental assumption, that faults are rare (remember that the fault rate is in the order of 10%) and homogeneous. The inferior performance of the typical semi-supervised method (µ-SVM) illustrates the well-known smoothness versus locality argument for deep learning [START_REF] Bengio | Scaling learning algorithms towards ai[END_REF][START_REF] Bengio | Representation learning: A review and new perspectives[END_REF]: the difficulty to model the highly varying decision surfaces produced by complex dependencies involving many factors.

As shown in the score distribution of Fig. 11, the proposed architecture of the CNN model separates anomalous layers significantly. This allows for great flexibility in choosing the working point for deployment in production in the CMS DQM. Depending on the cost of type 1 and type 2 errors for the detector operators the threshold can be set anywhere in [0.1, 0.9] score range. When using the CNN for the selection of good samples for training the regional algorithms, the working point is chosen not to favor specificity nor sensitivity, with a threshold equal to score 0.5.

The production algorithm targets a specific failure scenario of dead regions and produces a chamber-wise goodness assessment, without being capable of identifying a specific problematic layer in the chamber. For this reason we can not directly compare its performance with our local approach. For the sake of benchmarking our approach, we use our per-layer ground truth to label as bad any chamber with at least one problematic layer. We then ask the production algorithm if it indicates there is at least one faulty layer in a chamber. With this per-chamber label, we are able to estimate the specificity of the production algorithm to 91%, with a sensitivity of only 26%.

Another difference of our approach with respect to the production one is the performance with low statis- tics i.e. at the beginning of a run. As seen in Fig. 12, our CNN model gradually adds alarms until reaching stability. The production algorithm has the opposite behavior, generating a substantial fraction of false alarms in the early stages of the run.

Regional approach

To assess the performance of a given ensemble of channels we take as anomaly indicator the quantity:

k i = 1 j j ( ẋk i,j -ẍk i,j ) 2 ,
i.e., the MSE between the original sample given as input to the encoder ( ẋk i,j ) and the output of the decoder (ẍ k i,j ). With the objective of identifying the problematic region of the chamber, we exploit the granularity of the autoencoder information computing the MSE values for different set of channels. For example, we can compute the MSE for all the channels corresponding to a given read-out electronic board or, alternatively we can compute it per layer when tackling potential failures of the voltage distribution system.

We use this figure of merit on the sample with different voltage settings described in Section 6. Figure 13 shows good performance of all models, especially convolutional autoencoder. The distributions of the MSE for a well behaving and a problematic layer are shown in Fig. 14. The MSE distribution for layer 9 shows clear separation for chambers operated at nominal and lower voltages. For each i value for a given example, a quantitative assessment of the severity of a potential anomaly can be derived quoting the corresponding p-value of the Fig. 13 ROC and AUC of the different autoencoder models used in regional approach. The discriminator between good and anomalous samples is the in layer 9.

good example distribution. The separation is less pronounced for the working point at 3450 V being closer to the nominal setup. This reflects in the AUC values reported in Fig. 13.

The production algorithm is not sensitive to the type of faults described in this section since the hits in layer 9 are non-zero values. Thinking about deployment in the DQM infrastructure of the CMS experiment, the best result would be obtained when applying the local and regional models in a pipeline.

Global approach

Figure 15 shows an example of a low-dimensionality representation of the chamber data clustering depending on the chamber position in the detector. The global approach is then potentially capable to spot an unusual behavior of DT chambers taking into account the geographical constraints. Ultimately, this could pave the way to more flexible assessment by scoring per detector region.

When investigating the representations for a specific chamber across different runs (see Fig. 16), we notice that the representations tend to cluster depending on the number of problematic layers. Thanks to this fact, the cumulative distribution of the compressed representation could be used to highlight the occurrence of new anomalies or to associate an anomalous behavior to an already known problem. This application could assist experts in diagnosing transient and reoccurring issues.

Conclusions and Outlook

This paper shows how detector malfunctions can be identified with high accuracy by a set of automatic procedures, based on machine learning. We considered the specific case of the DT muon chambers of the CMS experiment. We developed a CNN-based classifier to spot local misbehaviors of the kind currently targeted by the existing CMS monitoring tools. We also showed that it is possible to extract more information from the map of electronic hits than the currently implemented statistical tests. In particular, we developed a strategy to spot regional problems across layers in a detector chamber, or globally, i.e., across chambers in the full muon detector. These algorithms, based on autoencoders, will offer a more robust anomaly detection strategy, not being defined as supervised classifiers of specific failure modes. This approach allows to localize the origin of a given anomaly, exploiting the granularity offered by the use of MSE of the decoded image as a quantification of the anomaly.

Currently, these algorithms have been integrated into the CMS online DQM infrastructure and they are being commissioned with the early data of the 2018 Run. The model could be further refined, e.g. integrating a mechanism of periodic retraining that would allow to repeat alarms for known problems, or to adapt to the long-term changes of the detector and beam conditions.

Since CNN is the basic ingredient in this study, and since many monitored quantities in typical high-energy physics experiments are based on 2D maps (e.g., detector occupancy, detector synchronization, etc.), the approach proposed in this paper could be extended beyond the presented use case. We hope that this case study could serve as a concrete showcase and could motivate further DQM automation using machine learning.

ing the computing resources to train our models and of CERN OpenLab for sponsoring A.S.'s internship at CERN, as part of the CERN OpenLab Summer student program. We thank Danilo Rezende for precious discussions and suggestions. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement n o 772369).

Fig. 1

 1 Fig. 1 Schematic view of the one DT chamber showing the position and orientation of the tubes. From [8].

Fig. 2

 2 Fig. 2 Magnified view of the CMS detector showing the wheel structure. The muon chambers are represented by the white volumes while the red volumes represent the iron return yoke.

Fig. 3

 3 Fig.3Numbering schema of the sector and stations of DT chambers in one wheel. From[START_REF]Calibration of the cms drift tube chambers and measurement of the drift velocity with cosmic rays[END_REF].

Fig. 4

 4 Fig. 4 Example of visualization of input data for three DT chambers. The data in (A) manifest the expected behavior in spite of having a dead channel in layer 1. The production algorithm regards this instance as non-problematic. The chamber shown in the plot in (B) instead shows regions of low occupancy across the 12 layers and should be classified as faulty. According to the run log, this effect was induced by a transient problem with the detector electronic. (C) suffers from a region in layer 1 with lower efficiency, which should be identified as anomalous. The production algorithm classify the chamber in (B) as anomalous. However it is not sensitive enough to flag the chamber in (C).

Fig. 5

 5 Fig. 5 Example of two kinds of input sample preprocessing. (A) acquired (raw) values, (B) standardizing each layer directly from raw values using linear interpolation. (C) smoothing the raw values data with median filter (D) standardizing each layer from smoothed data. In (C), the isolated lowoccupancy spot in layer 1, corresponding to a dead channel, is discarded.

Fig. 6

 6 Fig. 6 Architecture of the convolutional neural network model used to target the local strategy.

Fig. 7

 7 Fig. 7 Example of the impact on hit counting of different voltages applied to layer 9. (A) shows the occupancy map when operating the layer at 3200 V and (B) shows the effect of operating at 3450 V. Both examples should be regarded as anomalies. Since the values in both cases are not equal to zero, the algorithm considers those cases as nonproblematic.

Fig. 8

 8 Fig.8Convolutional (A) and simple, denoising, sparse (B) autoencoder models architecture used to target regional strategy.

Fig. 9

 9 Fig. 9 ROC curves for different models used in the local approach. The Area Under the Curve (AUC) is quoted to compare the performance.

Fig. 10

 10 Fig. 10 Loss function as a function of the number of epochs in the training of the CNN model used for the local approach. The two curves illustrate the behavor of the training and validation data sets.

Fig. 11

 11 Fig. 11 Distribution of scores in local approach for the CNN model.

Fig. 12

 12 Fig. 12 Stability of the proposed model and the production algorithm as a function of time (number of lumisections) for three different runs: 306777, 306793, 306794. The stability test is performed every 10 LSs. The CNN: Total and Production lines follow the total fraction of alarms generated by the methods. Instead CNN: Emerging reports the fraction of new alarms being generated by the CNN model w.r.t previous test point.

Fig. 14

 14 Fig. 14 MSE between reconstructed and input samples for layer 3 (A) and layer 9 (B) for 3 categories of data for convolutional autoencoder. Despite a problem in layer 9, all for layer 3 are comparable for all chambers.

Fig. 15

 15 Fig. 15 Compressed representation of the chamber-level data of the global model. The samples cluster according to position in the detector. Here depending on the station number.

Fig. 16

 16 Fig. 16 Compressed representation of the chamber-level data of the global model limited to only one chamber across different runs with respect to number of faulty layers (scale). The samples cluster according to similar behavior.

Acknowledgements

We thank the CMS collaboration for providing the data set used in this study. We are thankful to the members of the CMS Physics Performance and Data set project and the CMS DT Detector Performance Group for useful discussions, suggestions, and support. We acknowledge the support of the CMS CERN group for provid-


Xavier Mougeot

Reliability of usual assumptions in the calculation of β and ν spectra

Reliability of usual assumptions in the calculation of

β and ν spectra

I. INTRODUCTION

Beta emission properties are crucial features in nuclear data evaluations for building the decay schemes of radionuclides. Regarding β spectra, DDEP evaluations [1] as well as NNDC evaluations [START_REF] Bhat | Evaluated Nuclear Structure Data File (ENSDF)[END_REF] provide mean energies calculated by the LOGFT program [3]. The other code commonly used for nuclear data is the RADLST program [START_REF] Burrows | The Program RADLST[END_REF]. Both codes use the same assumptions, which are tested in the present work, and can calculate only allowed, first and second forbidden unique transitions. These assumptions are still used in the most recent calculations [START_REF] Huber | [END_REF]6].

Obviously, the mean energies depend directly on the shapes of β spectra. Furthermore, the energy distribution of the emission probability is needed in ionizing-radiation metrology [7], nuclear medicine for radiotherapy and dosimetry [8], as well as in the nuclear industry for calculations of residual power or post-irradiation fuel management [START_REF] Nichols | Summary Report of Consultants' Meeting on Beta-Decay and Decay Heat[END_REF]. The measurement of ν spectra from reactors, correlated with β spectra, is currently developing as a technique of monitoring nuclear facilities and would enable direct access to the amount of fissile material produced for military use [START_REF] Dos Anjos | Final Report: Focused Workshop on Antineutrino Detection for Safeguards Applications[END_REF]. The β decays are also known to play an important role in astrophysics, e.g., for the r process, in cosmochronology or for galatic neutrinos [START_REF] Klapdor | [END_REF]. In elementary particle physics, measured ν spectra could highlight possible physics beyond the standard model leading to a hypothetical fourth neutrino coupled only to the gravity field [12]. However, a recent study has shed light on the major influence of the uncertainties that can be associated with β-spectra calculations, casting doubt on the reality of the "reactor neutrino anomaly" [13].

Atomic effects have already been demonstrated in previous work having a substantial influence on the spectrum shape at low β energies, which correspond to high ν energies [14]. The calculations presented here are less refined. The aim of this work is to clearly point out the usual assumptions that are used in β calculations and to test them via a systematic comparison with recommended experimental shape factors. This study * xavier.mougeot@cea.fr is focused on simple and distinct transitions, disregarding measurements where multiple β emissions were not isolated or (neutrinoless) double β decays. The present calculations of the β spectra, performed using Behrens' formalism [15], are set out in detail in Sec. II. Regarding their importance for some research communities, the elements required to calculate the corresponding ν spectra are also given, but obviously the comparison with the measurements presented in Sec. III was performed only for β spectra.

II. CALCULATION OF β AND ν SPECTRA

Before going into the details of the calculation, some quantities used throughout this work have to be defined. The fine structure constant is denoted by α, while m e and m p stand for the electron and proton rest masses, respectively. The values of these constants are from the 2010 CODATA recommended values [START_REF] Mohr | [END_REF]. The ν particles are assumed to be massless in accordance with [15]. Indeed for this study, a nonzero mass less than 0.5 keV would affect very weakly the shape of β spectra in the endpoint region, where the emission probability is low. The daughter atomic number is Z and the mass number of the considered isobaric transition is A. The nuclear radius R of the daughter atom is given, according to [START_REF] Behrens | Numerical Tables for Beta-Decay and Electron Capture[END_REF], by R = 0.42587αA 1/3 , which corresponds to the usual formulation R = 1.2A 1/3 fm. Note that the present calculations are not very sensitive to this value. Other models, as the one given in [START_REF] Huber | [END_REF]6], were tested and did not lead to significant changes.

The particle parameters are the following. The kinetic energy of the β particle is denoted by E. The transition energy E 0 is the maximum energy of the β spectrum. The total energy of the β particle including its rest mass is W = 1 + E/m e . Similarly, the total energy corresponding to E 0 is

W 0 = 1 + E 0 /m e . The β-particle momentum is denoted by p = √ W 2 -1.
Finally, the ν-particle energy is q = W 0 -W , which is also its momentum because of the hypothesis of massless ν particles.

Consistent with Behrens' formalism, natural units = m e = c = 1 are used throughout this work. All numerical integrations were performed using a local three-point Lagrangian interpolation. Each spectrum was normalized to unity.

A. Spectrum shape

The shape of a β spectrum is proportional to (i) a phase space factor pW q 2 for the sharing of the momentum between the leptons, (ii) the so-called Fermi function F Beh. (Z,W ) which accounts for the static Coulomb corrections from the electromagnetic field of the nucleus, and (iii) a shape factor C(W ) which accounts for the remaining energy dependence that comes from the coupling between the nuclear structure and the leptons dynamic. Thus, following Behrens' formalism [15], it is given by

dN dW ∝ pW q 2 F Beh. (Z,W )C(W ). (1)
It should be highlighted that Behrens' formalism exhibits a certain symmetry for which the calculation of a β -decay is performed for a positive Z, while a negative Z is used for a β + decay. It is also noteworthy that in this formalism F Beh. (Z,W ) is defined differently from the usual Fermi function. Actually, it is the product of the usual Fermi function F 0 and a factor L 0 which accounts for the distortion of the electron wave functions by a non-point-like nuclear charge. F 0 L 0 can be determined directly from the electron wave functions (see Sec. II B). The theoretical shape factors C(W ) and the usual approximations tested in the present work (λ k = 1 and ξ approximation) are described in Sec. II C. The spectrum given in Eq. ( 1) is also corrected for the atomic screening effect (Sec. II D) and for radiative corrections (Sec. II E).

The present calculations do not take into account any weak magnetism correction because currently, there are neither theoretical nor experimental bases to establish a clear and trustworthy correction. A recent attempt proposes a first-order correction, but with an associated uncertainty much greater than the corrected effect [START_REF] Huber | [END_REF]6]. One can expect that a careful study of the experimental shape factors recommended in this work, especially for the allowed and first forbidden unique transitions, could lead to a more precise estimate of the weak magnetism effect. But such a study would require sophisticated calculations that account for the atomic and nuclear structures sufficiently well to clearly deconvolve them.

Because of the conservation of the leptonic charge and of the energy during the transition, there is a symmetry between β and ν spectra, the latter being the complement of the former. Indeed, the probability of a β particle to be created with energy W is the same as the probability of a ν particle to be created with energy q = W 0 -W [15]. The only awkward problem concerns the radiative corrections, which are different and lead to a small asymmetry (see Sec. II E).

B. Atomic potential and electron wave functions

Electron wave functions have to be calculated to determine some parameters (F 0 L 0 and λ k ) involved in the spectrum shape. Obviously, the structure of the Coulomb potential highly influences these wave functions. A more complicated potential, which was used to calculate precisely the atomic effects in allowed transitions, has already been described elsewhere [14]. Unfortunately, the required computation time for one spectrum is non-negligible and moreover, for most measurements referred to in this work the lowest energies recorded are too high to be sensitive to these atomic effects.

Therefore, the Coulomb potential of a uniformly charged sphere with a nuclear radius R was considered, and it was assumed to be scalar, static, and spherically symmetric. It exhibits a quadratic behavior inside the nucleus and is the same as for a pointlike charge outside,

V (r) = ⎧ ⎨ ⎩ -αZ 2R 3 -r R 2 0 r R, -αZ r R r ∞. (2)
The finite nuclear size effect is thus intrinsically included in the present calculations.

The relativistic behavior of the β particle within the Coulomb field of the nucleus is described by the Dirac equation. The electron wave function can be split into a radial part and an angular part (see [14] for more details). Thus, it is only necessary to calculate the electron radial wave functions f κ (r) (small component) and g κ (r) (large component). Here κ is the eigenvalue of the operator K = β( σ • L + 1) which appears by applying the theory of angular momentum to an electron in a Coulomb central field (β is the Dirac matrix, σ designates the Pauli matrices σ x,y,z , and L is the orbital angular momentum operator). Analytical solutions of the Dirac equation for f κ and g κ exist only for pure Coulomb potential, i.e., V (r) = -αZ/r [START_REF] Rose | Relativistic Electron Theory[END_REF]. The somewhat simple potential in Eq. (2) requires determining the radial wave functions numerically.

Bühring's method was followed, in which the electron radial wave functions are expressed locally as power series expansions, allowing the Dirac equation to be solved with recurrence relations [START_REF] Bühring | [END_REF]. Specific solutions have to be determined relative to the singular points r = 0 and r = ∞. Behrens and Bühring describe two methods in [15] for determining the global solution. The first one is called "the simple but slow method." The radial wave functions are evaluated (i) near the singular points with their specific solutions and (ii) step by step between these two points using the ordinary solution to (iii) reconnect each solution with the appropriate renormalization and phase shift. The so-called Coulomb amplitudes α κ are given, and defined, by the reconnection of the solutions near the origin with the ordinary solutions. These amplitudes are linked to the values of the radial wave functions at the nuclear surface. Similarly, the phase shifts are given by the reconnection of the asymptotic solutions with the ordinary solutions. This procedure is time consuming because it requires evaluating the wave functions in the whole space, without neglecting any spatial oscillation.

The second method, used in this work, is referred to as "the fast but complicated method." Actually, it is possible to avoid evaluating the wave functions in the whole space in the specific case of the potential given in Eq. ( 2). This method is too complicated to be set out in detail here, however, the basic idea is to choose the solutions relative to the origin as a basis for representing the solution outside the nucleus. The reader is urged to look at [15] for a complete explanation. Because of the complexity of this method, validation of the present calculations was attempted here. The unscreened Coulomb functions (Table II in [START_REF] Behrens | Numerical Tables for Beta-Decay and Electron Capture[END_REF]) were recalculated and perfectly reproduced up to the last decimal. Thus, the code implemented for this work can calculate these parameters at any energy, especially the Fermi function F 0 L 0 and the λ k parameters, avoiding any interpolation in these tables. Indeed in Behrens formalism, the Fermi function includes the finite nuclear size effect through the L 0 factor and is determined using the Coulomb amplitude of the electron wave functions,

F Beh. (Z,W ) = F 0 L 0 = α 2 -1 + α 2 1 2p 2 , ( 3
)
disregarding the structure of V (r) and the method used to calculate the wave functions. The λ k parameters are defined in Sec. II C.

C. Theoretical shape factors

The β decays are classified according to the spin change J = |J i -J f | and the parity change π i π f , where (J i ,π i ) and (J f ,π f ) denote the spins and parities of the initial and final nuclear states, respectively.

For allowed and forbidden unique transitions, Behrens and Bühring have demonstrated that in a first approximation, but with excellent precision, the nuclear current component can be factored out and only the dynamical factor from the coupling with the leptons remains [15]. Introducing L = 1 for J = 0 and L = J otherwise, the shape factor of an allowed (L = 1) or an (L -1) th forbidden unique transition is

C(W ) = (2L -1)! L k=1 λ k p 2(k-1) q 2(L-k) (2k -1)![2(L -k) + 1]! , (4)
with k = |κ| and λ k (Z,W ) is determined from the Coulomb amplitudes α κ defined in Sec. II B,

λ k = α 2 -k + α 2 k α 2 -1 + α 2 1 . ( 5
)
As described in Sec. II B, the procedure to calculate these λ k parameters is not straightforward. For this reason, a typical, but not justified, assumption in classical β-spectra calculations, even in the most recent ones [START_REF] Huber | [END_REF]6], is to set λ k ≡ 1. This assumption is referred to as the λ k = 1 approximation in this work.

However, calculation of β spectra for nonunique transitions is far more difficult than for unique transitions because the structures of the initial and final nuclear states have to be taken into account. It can be justified for some first forbidden nonunique transitions to calculate them as allowed if the ξ approximation is fulfilled [START_REF] Schopper | Weak Interactions and Nuclear Beta Decay[END_REF], i.e., when the Coulomb energy of the β particle emitted at the nuclear surface is much higher than the maximum energy of the transition,

2ξ = αZ R E 0 . ( 6
)
However, despite no theoretical justification, it is often generalized to all nonunique transitions, calculating a forbidden nonunique transition as a forbidden unique one having the same J . No other study was found that evaluates the applicability of the ξ approximation or that sets a limit on the ratio 2ξ/E 0 for which it would definitely provide good results.

D. Atomic screening effect

Screening is the only atomic effect taken into account in this work. Excitation and ionization of atomic electrons, respectively, called shake-up and shake-off effects, can arise from the sudden change of the nucleus charge during the β transition. These secondary effects can be confidently set apart because the probability of K-shell depletion per β decay is expected to be about 0.1% [START_REF] Cooper | [END_REF].

Studying the 63 Ni and 241 Pu decays, it was already proved in previous work that two main atomic effects have to be accounted for to accurately reproduce β spectra at low energy, the screening and the exchange effects [14]. A new and precise screening correction was defined at that time for allowed transitions that takes into account the spatial extension of both the potential V (r) and the wave functions, but its calculation is time consuming. In the present work, the common approach described below and originally set out 80 years ago by Rose [22] was preferred to test its influence. This correction is largely found in the literature and still used in modern codes.

In this approach, the screening effect is approximately estimated by subtracting a screened potential of a Thomas-Fermi type that depends only on Z. For β -transitions, the potential was taken from [START_REF] Dillman | EDISTR -A Computer Program to Obtain a Nuclear Decay Data Base for Radiation Dosimetry[END_REF] V 0 (β -) = (-9.45 × 10 -9 )|Z| 3 + (3.014 × 10 -6 )|Z| 2 +(1.881 × 10 -4 )|Z| -(5.116 × 10 -4 ).

(

) 7
For β + transitions, the potential is modified in the following way [START_REF] Gove | [END_REF]:

V 0 (β + ) = V 0 (β -) exp(a/p + b/p 2 ), (8)
with the parameters,

a = (1.11 × 10 -7 )|Z| 3 -(1.01 × 10 -5 )|Z| 2 -(2.38 × 10 -3 )|Z| + 0.102, (9)
and

b = (-2.42 × 10 -8 )|Z| 3 + (3.83 × 10 -6 )|Z| 2 + (3.60 × 10 -5 )|Z| -0.0156. ( 10
)
The method for applying this correction is well described in [START_REF] Gove | [END_REF] and consists of evaluating the probability of β - emission at a given energy W as the same probability as at the energy W = W -V 0 . This can be summarized by a simple three-part rule, where the screening is applied everywhere in the spectrum shape in Eq. ( 1) and in the subsequent factors in Eqs. (3)-( 5) by the following: replacing W by W ; replacing p by the corresponding momentum p = √ W 2 -1; keeping the ν energy q as it is. Obviously, this correction cannot be applied for W -V 0 < 1 because p would become imaginary. This proviso induces a nonphysical discontinuity in the spectrum at the energy E sc = m e V 0 , where the screening correction can be switched on. Conversely for β + transitions, the energy W is replaced by W = W + V 0 , the momentum p by the corresponding p , and the correction can be applied across the entire energy range of the spectrum.

E. Radiative corrections

Up to now, only static Coulomb corrections have been considered and atomic effects have been taken into account only through the screening correction, excluding all other atomic processes (see Sec. II D). The radiative corrections are nonstatic Coulomb corrections and they can have an important effect on β-spectra shapes for high E 0 . Formally, they should be determined using Feynman diagrams, which is an approach far outside the present work. However, they can be split into "inner" and "outer" radiative corrections and only the latter influences the energy dependence of the β spectrum [25]. Thus, the "inner" radiative corrections have been left out in present calculations.

Outer radiative corrections account for the energy loss of the β particle in the electromagnetic field of the nucleus by means of virtual photons and internal bremsstrahlung. In the present work, they are calculated as described in [15,25,26]. The spectrum is corrected by a factor [1 + δ R (W,Z)] defined as follows:

δ R (W,Z) = δ 1 (W ) + δ 2 (Z) + δ 3 (Z) + δ 4 (Z). ( 11
)
The energy-dependent term is of order α,

δ 1 (W ) = α 2π g(W,q), ( 12
)
-R 2 ), ( 1
-R 2 ) λ k =1 , Ē
and Ēλ k =1 are in %.

Experimental shape factor

Parent Mode

E 0 ×(q 2 + λ 2 p 2 ) except for 204 Tl Ref. Range %E 0 Ēsf Ēcalc (1 -R 2 ) Ē Ēλ k =1 (1 -R 2 ) λ k =1 Ēλ k =1
First forbidden unique where

g(W,q) = 3 ln m p m e - 3 4 + 4 β L 2β 1 + β + 4 tanh -1 β β -1 q 3W - 3 2 + ln (2q) + tanh -1 β β (1 + β 2 ) q 2 3W 2 -4tanh -1 β , ( 13
)
with β = p/W and L(x) is the Spence function, which is the opposite of the dilogarithm function,

L(x) = x 0 ln(1 -t) t dt. ( 14
)
The higher order terms in (αZ) are to a first approximation energy independent,

δ 2 (Z) = 1.1 |Z|α 2 m p m e , ( 15
)
δ 3 (Z) = Z 2 α 3 π 3 ln 2 - 3 2 + π 2 3 m p m e , ( 16
)
δ 4 (Z) = |Z|α 3 2π m p m e . ( 17
)
These radiative corrections have recently been calculated for ν particles and it was proved that the internal bremsstrahlung contribution is different, the rationale being that nonstatic Coulomb processes are secondary processes [27]. Thus, in present calculations of ν spectra, only the following result from [27] is accounted for, which modifies the spectrum by a factor [1 + δ ν (q)] defined as follows:

δ ν (q) = α 2π h(W ), (18)
with

h(W ) = 3 ln m p m e + 23 4 + 8 β L 2β 1 + β + 8 tanh -1 β β -1 ln(2Wβ) + 4 tanh -1 β β 7 + 3β 2 8 -2tanh -1 β , ( 19
)
where W is still the β-particle energy that corresponds to the ν energy q, and the quantity β = p/W still remains the same. This difference in the radiative corrections is the only one in the present work that breaks the symmetry between a given β spectrum and its corresponding ν spectrum.

III. SYSTEMATIC COMPARISON WITH EXPERIMENTAL SHAPE FACTORS

A database of 130 β transitions was created from the experimental shape factors available in the literature. These shape factors come mainly from [START_REF] Behrens | Shapes of beta Spectra[END_REF] but with a different choice of the recommended publication in some cases. The literature since [START_REF] Behrens | Shapes of beta Spectra[END_REF] was also reviewed and new shape factors have been considered. The vast majority (72%) were established from measurements using magnetic spectrometers. It is noteworthy that for some transitions, significant disagreement exists between different publications. Each transition energy was updated with the latest DDEP evaluation [START_REF]Recommended Atomic and Nuclear Data Evaluated by the DDEP Collaboration[END_REF], or with the latest NNDC evaluation [START_REF] Bhat | Evaluated Nuclear Structure Data File (ENSDF)[END_REF] if the parent radionuclide was not yet evaluated by the DDEP collaboration.

Calculations have been performed as described in Sec. II to analyze the adequacy of these experimental shape factors. Each spectrum was calculated across the full energy range of the transition. The analysis was handled using parameters defined in Sec. III A and the results are extensively presented in Sec. III B, grouped by the transition type. The calculations have been performed with and without the screening correction presented in Sec. II D. The influence of this correction on the results is weak, especially when E 0 is high, but no systematic trend was found. Therefore, only results with screening accounted for are given in the present publication.

A. Analysis

For each β transition, the experimental shape factor was not necessarily measured in the full energy range. It is an important feature of the measurement which clearly has a significant influence on the present analysis. Thus, energy range of the measurement and proportion (%E 0 ) of this range compared to the maximum available energy E 0 are provided for each considered transition. This simple criterion highlights that in many cases, the β spectra are not so well measured: Forty-four

Ēsf Ēcalc (1 -R 2 ) Ē Ēλ k =1 (1 -R 2 ) λ k =1 Ēλ k =1
Second forbidden unique a See [14] for more precise results.

transitions over a total of 130 exhibit (%E 0 ) < 50%, namely about one-third of the transitions. The mean energy of each calculated β spectrum, from an experimental shape factor or according to one or another assumption, was determined as follows:

Ē = m e W 0 1 (W -1)N (W )dW W 0 1 N (W )dW . ( 20
)
The difference Ē with the mean energy Ēsf determined from the experimental shape factor is also calculated.

To evaluate the global quality of present calculated spectra, a goodness-of-fit parameter was estimated. Considering a set of n measurements {y i } and a model that leads to n predictions {y th i }, the quality of the fit of this model with respect to the data can be given by

R 2 = 1 - var( êi ) var(y i ) . ( 21
)
This is a general definition which is valid for n 1, without any restriction on the type of model. This quantity falls in the range 0 R 2 1, where R 2 = 1 corresponds to a perfect prediction. The modeling error is êi = y iy th i . Nonbiased estimators of the variances are

var(y i ) = n i=1 (y i -y) 2 n -1 , var( êi ) = n i=1 ( êi ) 2 n -N p -1 , (22)
with y the simple mean of the measurements, and N p the number of model parameters. For β spectra, E 0 is clearly one parameter. An additional one is also needed to take into account that the shape factors had been determined from measured spectra, which somewhat hides a normalization. Thus, N p = 2 was considered to determine R 2 . Therefore, the parameter (1 -R 2 ) quantifies the general disagreement between the measurements and the predictions. For the consistency of the analysis, the (1 -R 2 ) values have only been calculated in the energy range of the measurements.

For forbidden unique transitions, the mean energy Ēcalc refers to the calculated spectrum for which the λ k = 1 approximation was not considered. This approximation was used to calculate the mean energy Ēλ k =1 and all other parameters having the same label. For forbidden nonunique transitions, the ratio 2ξ/E 0 is also given because it should determine, in principle, the validity of the ξ approximation. The mean energy Ēcalc and all other parameters then refer to the calculated spectrum according to the ξ approximation, as defined in Sec. II C.

In the following, a β spectrum is considered well known in the energy range of the measurement if (1 -R 2 ) < 1% and | Ē| < 1%, which will be referred to as the 1% criterion hereafter. If (1 -R 2 ) is low but | Ē| is high, it means that the calculation is correct only inside this energy range and, because the mean energy is determined using the full energy range of the transition, then the agreement is significantly worse outside.

B. Results

The results of this systematic comparison with experimental shape factors are presented first for the allowed transitions in Sec. III B 1. Next, they are shown for the forbidden unique transitions in Sec. III B 2 and the reliability of the λ k = 1 approximation is reviewed. Finally, the results for the forbidden nonunique transitions are given in Sec. III B 3 and a critical point of view on the ξ approximation is provided.

Allowed transitions

Results for the allowed transitions are presented in Table I.

The disagreement parameter (1 -R 2 ) is found to be less than 0.45% in all cases. With | Ē| < 1%, 30 of the 36 tabulated spectra are reproduced well. It should be noted that better results were found for the 63 Ni decay with more precise calculation of the atomic effects [14]. If this criterion is Only the 14 C decay remains for which | Ē| = 2.51%, which does not fail so badly. This allowed transition exhibits an anomalously long half-life, which clearly suggests a nuclear structure effect which would also affect the spectral shape [START_REF] Towner | [END_REF].

Forbidden unique transitions

In Table II are presented the results for the first forbidden unique transitions, and in Table III for the second and third forbidden unique ones.

Using the 1% criterion defined in Sec. III A, 21 of the 25 first forbidden unique spectra are reproduced well. Extending the disagreement parameter to (1 -R 2 ) < 2% and the mean energy parameter to | Ē| < 2%, 24 transitions can be seen as well known. The only one for which the calculation really fails is the 166 Ho decay. Despite the measured energy covers 72% of E 0 and an emission probability of 50%, the quality of this measurement is perhaps questionable because the authors did not provide much detail of possible source impurities or the data analysis [62]. A new precise measurement would be beneficial.

It is difficult to be conclusive regarding the second and third forbidden unique transitions because of the very small number of measurements. With the 1% criterion, three of the four second forbidden unique spectra are reproduced well. The last one is the 138 La decay which was very precisely measured in recent work down to 2.5 keV and for which atomic effects are expected to play an important role [63]. The third forbidden unique transition from the 40 K decay exhibits a disagreement parameter (1 -R 2 ) = 1.21% and | Ē| = 0.122%. As it is the only measured transition of this type, it cannot be established if the present calculations are correct for this type of transition in general.

From this analysis, it was found that for all the 30 forbidden unique transitions, results from the λ k = 1 assumption are systematically worse than the full calculation with the correct λ k except for the first forbidden unique decay of 41 Ar and the second forbidden unique decays of 22 Na and 10 Be. In all three transitions the results pass the 1% criterion. Thus, these isolated cases cannot be rationally considered as significant regarding the relevance of the λ k = 1 approximation.

Forbidden nonunique transitions

Results are shown in Table IV for the first forbidden nonunique decays, and in Table V for the second, third, and fourth forbidden nonunique ones. Analysis parameters (1 -R 2 ) and | Ē| are presented in Fig. 1 as a function of 2ξ/E 0 .

Applying the 1% criterion, 26 of the 53 first forbidden nonunique transitions are reproduced well, namely 49%. For the 241 Pu decay, better results were found with a more precise calculation of the atomic effects [14]. None of the second and third forbidden nonunique transitions are well reproduced, but the fourth one from the 113 Cd decay passes the 1% criterion. However, no clear statement about the validity of the ξ approximation for this transition type can be made from this single measurement. It is difficult to be conclusive regarding the ξ approximation for extricating a correct criterion. Indeed, one can see in Fig. 1 and in Tables IV and V that the criterion 2ξ/E 0 1 is not sufficient to predict the correctness of the calculations. A conservative criterion could be 2ξ/E 0 > 100, which is confirmed by the first forbidden nonunique decay of 151 Sm, not present in the tables because no shape factor was yet extracted. However, this transition, with E 0 = 76.6 keV and 2ξ/E 0 = 204, was recently measured and proved to exhibit an allowed shape [83]. Below 2ξ/E 0 = 100, no systematic trend can be observed.

IV. CONCLUSION

Setting out in detail the β and ν calculations used in this work, the usual assumptions have been clearly pointed out. A database of measured shape factors from distinct β emissions was created. The usual assumptions were tested by comparison with the measurements and focus was made on the λ k = 1 approximation for the forbidden unique transitions and the ξ approximation for the forbidden nonunique ones. The main conclusions that can be drawn from this work are the following:

(1) For allowed transitions, except in the specific case of the 14 C decay, excellent agreement was found between theory and experiment. (2) The λ k = 1 approximation was proved to be systematically worse than the full calculation that includes the energy dependence of the λ k parameters. This assumption should be avoided in precise β and ν spectra evaluations.

(3) For the first time, the relevance of the ξ approximation was thoroughly studied for many transitions. About half of the tabulated first forbidden nonunique transitions were reproduced well, but significant disagreement was found in all other cases. No specific tendency arose from the analysis of the present results other than highlighting that the condition 2ξ/E 0 1 is not sufficient in itself. A conservative limit of the 2ξ/E 0 ratio of 100 was deduced, but has to be confirmed experimentally.

This study demonstrates also that the mean energy of a β spectrum given in the nuclear databases is definitely erroneous for forbidden transitions. And even for an allowed transition, this mean energy is not trustworthy in the case of low E 0 because of the major influence of atomic effects. A full calculation is thus recommended for each data evaluation.

However, improving nuclear databases, e.g., with a recommended spectrum or an experimental shape factor, will not be helpful if this information is not systematically taken into account. Presently, most research communities and the nuclear power industry use their own codes to calculate the mean energies from the transition parameters given in the databases. These codes are based on the same assumptions which have been tested in the present work and thus an effort should be made to upgrade them. Furthermore, in the context of reactor ν oscillations, a strong impact can be expected on the method used to convert measured total fission β spectra to total ν spectra. This impact should, at least, be estimated with improved calculations. Unfortunately, a typographical error is present in Eq. ( 13) of Sec. II E for the radiative corrections. Originally, the third line of Eq. ( 13) which defines the quantity g(W,q) was given as

• • • + tanh -1 β β (1 + β 2 ) q 2 3W 2 -4tanh -1 β ,
but should be

• • • + tanh -1 β β 2(1 + β 2 ) + q 2 6W 2 -4tanh -1 β ;
see Refs. [15,25]. Note that the correct formula was used throughout the work, and therefore all results given previously remain valid.

FIG. 1 .

 1 FIG. 1. (Color online) Distribution according to 2ξ/E 0 of the (1 -R 2 ) and | Ē| parameters for the measured shape factors of forbidden nonunique transitions.

Erratum:

  Reliability of usual assumptions in the calculation of β and ν spectra [Phys. Rev. C 91, 055504 (2015)]

TABLE I .

 I Experimental shape factors of allowed transitions. Energy range, Ēsf and Ēcalc are in keV; %E 0 , (1 -R 2 ) and Ē are in %.

	Parent Mode	E 0	Experimental shape factor	Ref.	Range	%E 0	Ēsf	Ēcalc (1 -R 2 )	Ē
	Allowed
	3 H	β -	18.591	1	[31]	11-18	37.7	5.696	5.695 2e-05 -0.009
	6 He	β -	3507.8	1 -0.2/W	[32]	300-3200 82.7 1585.30 1561.98	0.430 -1.47
	12 B	β -13369.3	1 + 0.0028W	[33,34] 7000-12 000 37.4 6479.16 6400.33	0.405 -1.22
	11 C	β +	960.5	1 -0.0074W	[35]	266-892 65.2 385.33 385.12	2e-04 -0.054
	14 C	β -	156.476	1 + 0.00124m e q	[36]	20-155	86.3	48.233 49.445 0.190	2.51
	12 N	β + 16316.1	1 -0.0027W	[33,34] 6000-13 000 42.9 7872.06 7870.28	0.004 -0.023
	13 N	β +	1198.45	1 + 0.0014W	[37]	100-1000 75.1 492.20 490.88	0.016 -0.267
	18 F	β +	633.9	1 + 0.0034W	[38]	70-600	83.6 250.50 249.71	0.006 -0.313
	20 F	β -	5390.86	1 + 0.006W -0.004/W -0.00009W 2	[39]	600-5380 88.7 2490.96 2468.62	0.120 -0.897
	22 Na	β +	546.44	1 -0.005W	[40]	70-520	82.4 216.32 215.68	0.004 -0.293
	24 Na	β -	1392.72	1 -0.011W	[41]	100-1390 92.6 552.78 553.87	0.004	0.196
	27 Si	β +	3789.83	1 -0.002W	[42] 1530-3800 59.9 1718.07 1713.32	0.006 -0.276
	31 Si	β -	1491.5	1 -0.003W	[43]	200-1300 73.8 594.93 594.64	0.002 -0.049
	32 P	β -	1710.66	1 -0.019W	[43]	200-1635 83.9 689.48 693.42	0.035	0.571
	34 Cl	β +	4469.64	1 -0.0003W	[44] 2500-4250 39.2 2052.22 2043.51	0.036 -0.424
	34m Cl	β +	2488.43	1 -0.0097W	[44] 1380-2300 37.0 1094.47 1095.36	0.003	0.081
	41 Ar	β -	1197.96	1 -0.017W	[45]	80-1140 88.5 456.47 458.34	0.018	0.409
	47 Ca	β -	1992	1 -0.11W -0.33/W + 0.009W 2	[46]	710-1899 59.7 817.07 816.96	0.046 -0.014
	49 Sc	β -	2006	1 -0.008W	[47]	150-1900 87.2 820.70 822.14	0.005	0.175
	56 Mn	β -	2848.7	1 + 0.29/W	[48]	970-2715 61.3 1195.02 1213.53	0.196	1.55
	56 Co	β +	1458.9	1 + 0.3/W	[49]	470-1330 59.0 622.98 629.74	0.135	1.09
	58 Co	β +	475.1	1 + 0.3/W	[50]	150-450 63.1 200.96 201.07	0.026	0.053
	60 Co	β -	317.32	1	[51]	75-280	64.6	95.52	95.50	1e-04 -0.030
	63 Ni a	β -	66.98	1 + 1213.49W -232.009/W	[52]	1-67	98.5	17.21	17.41	0.123	1.12
				-1517.16W 2 + 532.69W 3
	66 Ga	β +	4153	(1 + 0.03282W + 0.354/W -0.000419W 2 ) [53] 1550-4050 60.2 1882.36 1895.54	0.122	0.700
				×(1 + 0.4/W )
	68 Ga	β +	1899.1	1 -0.01 • W	[54]	870-1845 51.3 834.88 833.69	0.001 -0.143
	89 Zr	β +	902	1 -0.39W + 0.09W 2	[55]	100-840 82.0 392.22 395.03	0.271	0.718
	110m Ag β -	529.9	1 -0.01W	[56]	110-500 73.6 164.81 165.05	0.002	0.147
	115 Cd	β -	584.5	1 -0.0007W	[57]	250-570 54.8 184.70 184.74	3e-04	0.023
	115 Cd	β -	620.1	1 -0.0009W	[57]	250-620 59.7 197.89 197.96	4e-04	0.035
	114 In	β -	1988.6	1 -0.0015W	[58]	210-1930 86.5 777.17 777.38	0.001	0.027
	130 I	β -	1005	1 + 0.04W	[59]	750-980 22.9 349.88 346.63	0.311 -0.928
	130 I	β -	587	1 + 0.04W	[59]	170-460 49.4 185.35 184.23	0.030 -0.606
	131 I	β -	606.3	1 + 0.02W	[60]	380-590 34.6 191.90 191.31	0.017 -0.307
	134 Cs	β -	415.64	1	[61]	90-415	78.2 123.44 123.36	4e-04 -0.064
	134 Cs	β -	658.39	1	[61]	415-658 36.9 209.93 210.00	4e-04	0.033

a See

[14]

for more precise results.

TABLE II .

 II Experimental shape factors of first forbidden unique transitions. Energy range, Ēsf

	0 , (1
	and Ēλ k =1 are in keV; %E
	, Ēcalc

TABLE III .

 III Experimental shape factors of second and third forbidden unique transitions. Energy range, Ēsf , Ēcalc , and Ēλ k =1 are in keV; %E 0 , (1 -R 2 ), (1 -R 2 ) λ k =1 , Ē, and Ēλ k =1 are in %.

	Parent Mode	E 0	Experimental shape factor Ref.	Range	%E 0

TABLE IV .

 IV Experimental shape factors of first forbidden nonunique transitions. Energy range, Ēsf , and Ēcalc are in keV; %E 0 , (1 -R 2 ), and Ē are in %.

	Parent	Mode	E 0	Experimental shape factor	Ref.	Range	%E 0 2ξ/E 0	Ēsf	Ēcalc	(1 -R 2 )	Ē
	First forbidden nonunique
	42 K	β -	2000.8	1 -0.11W + 0.7/W + 0.012W 2	[84]	200-1850 82.5	3.45	769.42	822.39	5.27	6.88
	72 Ga	β -	3163.1	q 2 + λ 2 p 2 + 2.4	[ 85]	2540-3030 15.5	2.92	997.77 1353.51	4199	35.7
	72 Ga	β -	2533.2	q 2 + 0.95p 2	[86]	1970-2320 13.8	3.64 1066.96 1055.07	51.3	-1.11
	74 As	β +	944.6	1 -0.017W	[66]	100-880	82.6	9.68	409.07	407.32	0.023 -0.429
	74 As	β -	718	1 + 0.013W	[66]	100-600	69.6	13.53	243.00	242.43	0.004 -0.235
	76 As	β -	2402.9	1 -0.007W	[67]	1730-2270 22.5	4.01	988.95	990.94	0.006	0.202
	84 Rb	β +	783	1 -0.073W	[66]	100-665	72.2	12.60	340.13	340.75	0.119	0.183
	86 Rb	β -	699.2	1 -0.64W -0.46/W + 0.14W 2	[87]	140-615	67.9	14.77	234.61	232.81	0.711 -0.766
	91 Y	β -	339.5	1	[88]	90-290	58.9	31.43	100.27	100.21	3e-04	-0.068
	99 Mo	β -	848.1	1 + 0.049W -0.0047/W + 0.034W 2	[89]	470-780	36.6	13.15	298.48	289.36	1.94	-3.05
	99 Mo	β -	1214.5	1 -0.01W	[89]	450-1150 57.6	9.18	440.80	442.21	0.009	0.320
	111 Ag	β -	694.7	1 -0.01W	[90]	420-640	31.7	17.25	226.07	226.58	0.012	0.224
	111 Ag	β -	1036.8	1 -0.016W	[90]	700-990	28.0	11.56	361.72	363.48	0.051	0.487
	115m Cd	β -	339.1	1 + 0.011W	[73]	280-296	4.7	35.65	98.92	98.73	0.090 -0.201
	115m Cd	β -	1629.7	1 -0.78W -17.2/W	[91]	730-1580 52.2	7.42	533.21	617.82	23.5	15.9
	122 Sb	β -	1419.6	1 -0.012W + 0.003W 2	[92]	360-1380 71.9	8.86	522.73	522.76	0.003	0.004
	124 Sb	β -	2301.6	q 2 + λ 2 p 2 + 7.3	[ 93]	1680-2250 24.8	5.44	757.57	916.72	200	21.0
	124 I	β +	1534.9	1 -0.046W	[66]	820-1350 34.5	8.15	681.87	685.25	0.156	0.496
	139 Ba	β -	2151	1 -0.219W + 0.147/W + 0.02W 2	[71]	300-2060 81.8	6.14	752.68	840.60	7.54	11.7
	140 La	β -	2165.7	q 2 + λ 2 p 2 + 2.3	[ 94]	1680-2100 19.4	6.19	631.24	845.68	1470	34.0
	141 Ce	β -	435	1 -0.28W	[95]	170-420	57.5	31.27	122.55	129.51	2.28	5.68
	141 Ce	β -	580.4	1 -0.28W	[95]	170-550	65.5	23.43	167.10	180.53	2.63	8.04
	144 Ce	β -	318.6	1 -0.342W	[96]	174-300	39.6	42.39	86.33	91.07	5.75	5.48
	142 Pr	β -	586.6	1	[97]	120-510	66.5	23.52	182.48	182.52	5e-04	0.019
	143 Pr	β -	934.1	1 -0.018W	[98]	125-890	81.9	14.74	313.01	314.73	0.017	0.551
	144 Pr	β -	2997.4 1 + 0.0376W -0.118/W -0.0077W 2	[96]	350-2860 83.7	4.58 1219.55 1220.19	0.086	0.053
	147 Nd	β -	364.7	1 -0.2W	[99]	150-350	54.8	38.03	102.95	105.79	0.652	2.77
	147 Nd	β -	804.6	1 -0.17W	[100]	375-775	49.7	17.24	250.07	263.63	3.14	5.42
	146 Pm	β -	795	1 -0.394W -0.275/W + 0.044W 2	[101]	200-780	73.0	17.77	235.52	259.55	3.34	10.2
	150m Eu	β -	1015	1	[102]	80-970	87.7	14.24	344.42	344.77	0.002	0.103
	152 Eu	β -	1474.5	q 2 + 0.79p 2 + 5	[ 86]	1100-1450 23.7	9.76	532.81	534.76	4.40	0.367
	154 Eu	β -	1845.3	q 2 + 0.807p 2 + 12	[103]	1000-1770 41.7	7.76	689.89	694.75	0.587	0.705
	166 Ho	β -	1854.5 1 -0.87W -1/W + 0.22W 2 -0.02W 3	[96]	410-1790 74.4	8.01	672.36	693.20	0.396	3.10
	170 Tm	β -	883.7	1 + 0.125(q 2 + p 2 )	[ 104]	290-844	62.7	17.16	297.13	290.10	0.528 -2.37
	170 Tm	β -	968	1 + 0.103(q 2 + p 2 )	[ 104]	290-930	66.1	15.66	329.39	322.77	0.432 -2.01
	172 Tm	β -	1801	q 2 + λ 2 p 2	[75]	638-1740 61.2	8.39	440.64	667.28 108	51.4
	175 Yb	β -	356.3	1 -1.36W + 0.467W 2	[105]	220-320	28.1	42.75	61.21	102.12	10846	66.8
	176m Lu	β -	1224.8	1 -0.064W	[106]	250-1175 75.5	12.59	413.29	423.76	0.276	2.53
	176m Lu	β -	1313.2	1 -0.15W + 0.0167W 2	[106]	250-1280 78.4	11.74	442.44	459.79	0.631	3.92
	185 W	β -	432.5	1	[107]	90-370	64.7	36.52	126.59	126.48	0.002 -0.084
	187 W	β -	1310.9	q 2 + λ 2 p 2 + 23	[108]	720-1260 41.2	12.00	436.68	456.10	2.81	4.45
	186 Re	β -	932.3	1 + 0.038(q 2 + p 2 )	[ 104]	200-860	70.8	17.14	308.22	305.72	0.075 -0.813
	186 Re	β -	1069.5	1 + 0.034(q 2 + p 2 )	[ 104]	200-970	72.0	14.94	361.45	358.87	0.095 -0.715
	188 Re	β -	1965.3	1 + 0.031(q 2 + p 2 )	[ 104]	210-1840 82.9	8.10	731.69	728.54	0.419 -0.431
	188 Re	β -	2120.4	1 + 0.039(q 2 + p 2 )	[ 104]	210-1970 83.0	7.51	798.70	795.09	0.746 -0.453
	194 Ir	β -	1918.5	1 -0.01W	[109]	1300-1800 26.1	8.43	701.04	705.24	0.075	0.600
	198 Au	β -	961	1 -0.057W	[110]	100-956	89.1	17.14	308.77	314.44	0.145	1.84
	199 Au	β -	452	1 -0.24W	[95]	300-430	28.8	36.38	126.38	132.23	7.03	4.63
	206 Tl	β -	1532.4	1 -0.02W	[111,112] 195-1470 83.2	10.87	533.37	538.69	0.057	0.997
	207 Tl	β -	1418	1 + 0.024W	[113]	130-1290 81.8	11.73	495.65	492.32	0.033 -0.673
	209 Pb	β -	644	1 + 0.008W	[114]	200-613	64.1	26.06	197.17	197.00	5e-05	-0.087
	210 Bi	β -	1161.2	1 -0.47W + 0.065W 2	[115]	60-600	46.5	14.60	317.29	388.66	54.4	22.5
	241 Pu a	β -	20.8	1 -1.9582W + 0.96078W 2	[116]	1-21	96.2 880.64	5.10	5.24	0.291	2.75

TABLE V .

 V Experimental shape factors of second, third, and fourth forbidden nonunique transitions. Energy range, Ēsf , and Ēcalc are in keV; %E 0 , (1 -R 2 ), and Ē are in %.

	extended to | Ē| < 2%, 35 transitions can be considered as
	well known.
	Parent	Mode	E 0	Experimental shape factor	Ref.	Range	%E 0	2ξ/E 0	Ēsf	Ēcalc	(1 -R 2 )	Ē
	Second forbidden nonunique
	36 Cl	β -	709.53	1 -0.282715W -0.045988/W	[117]	6-710	99.2	9.22	314.09	278.40	25.0	-11.4
				-0.491739W 2 + 0.438731W 3
	46 Sc	β -	1477.2	q 2 + 1.67p 2	[81]	1170-1400	15.6	4.99	657.73	603.04	19.1	-8.31
	59 Fe	β -	1565	q 2 + 0.59p 2	[118]	470-1550	69.0	5.32	582.92	633.75	5.54	8.72
	58 Co	β +	1285.9	p 2	[119]	750-1150	31.1	6.27	708.71	580.91	114	-18.0
	94 Nb	β -	471.5	q 2	[120]	160-450	61.5	23.51	94.38	165.74	227	75.6
	99 Tc	β -	293.8	q 2 + 0.54p 2	[121]	55-250	66.4	38.84	95.28	101.03	2.72	6.04
	129 I	β -	151.2	3.16q 2 + p 2	[122]	55-145	59.5	84.81	45.98	48.62	2.88	5.75
	135 Cs	β -	268.7	q 2 + 0.10p 2	[123]	50-195	54.0	48.75	61.57	88.44	50.6	43.6
	137 Cs	β -	1175.63	see reference	[124]	710-1125	35.3	11.09	284.90	422.87	1834	48.4
	Third forbidden nonunique
	87 Rb	β -	283.3	q 4 + 0.305q 2 p 2 + 0.011p 4	[125]	65-180	40.6	36.32	56.65	114.99	262	103
	Fourth forbidden nonunique
	113 Cd	β -	322	3.034q 6 + 1.499 × 7q 4 p 2	[126]	30-320	90.1	37.77	139.83	139.30	0.076	-0.376
				+1.016 × 7q 2 p 4 + p 6

ACKNOWLEDGMENTS

This study was produced under the auspices and scientific support of the Decay Data Evaluation Project (DDEP) international collaboration.