A new software implementation of the Oslo method with rigorous statistical uncertainty propagation

TL;DR

This paper introduces OMpy, a new software implementation of the Oslo method that rigorously propagates uncertainties using Monte Carlo techniques, improving the analysis of nuclear level densities and gamma-decay strengths.

Contribution

The paper presents a novel software tool that enhances the Oslo method with comprehensive uncertainty propagation and simultaneous normalization of key nuclear parameters.

Findings

Successful verification with synthetic spectra

Improved uncertainty handling in nuclear data analysis

Comparison shows consistency with previous software

Abstract

The Oslo method comprises a set of analysis techniques designed to extract nuclear level density and average $\gamma$-decay strength function from a set of excitation-energy tagged $\gamma$-ray spectra. Here we present a new software implementation of the entire Oslo method, called OMpy. We provide a summary of the theoretical basis and derive the essential equations used in the Oslo method. In addition to the functionality of the original analysis code, the new implementation includes novel components such as a rigorous method to propagate uncertainties throughout all steps of the Oslo method using a Monte Carlo approach. The resulting level density and $\gamma$-ray strength function have to be normalized to auxiliary data. The normalization is performed simultaneously for both quantities, thus preserving all correlations. The software is verified by the analysis of a synthetic spectrum and compared to the results of the previous implementation, the oslo-method-software.

Figures

Equations

$$P(E_x,E_{\gamma })\propto \rho (E_f=E_x-E_{\gamma })\mathcal{T}(E_{\gamma }),$$

$$\mathcal{T}(E_{\gamma })=2\pi E_{\gamma }^{3}f(E_{\gamma }).$$

$$\rho (E_f=E_x-E_{\gamma }),f(E_{\gamma })\stackrel{G}{\to }A{e}^{\alpha E_f}\rho (E_f),B{e}^{\alpha E_{\gamma }}f(E_{\gamma }).$$

$${\mathcal{P}}_{\lambda }=p(k|\lambda )=\frac{{\lambda }^k{e}^{-\lambda }}{k!}$$

$$⟨{\vec{r}}_i⟩=\frac{1}{N_{\mathrm{ens}}}\sum _{l=1}^{N_{\mathrm{ens}}}{r}_i^{(m)},$$

$${\sigma }_i$$

$$P_{\mathrm{fit}}(E_x,E_{\gamma })=C_{E_x}\rho (E_x-E_{\gamma })\mathcal{T}(E_{\gamma }),$$

$${\chi }^2=\sum _{E_x,E_{\gamma }}{\left(\frac{P_{\exp }(E_x,E_{\gamma })-P_{\mathrm{fit}}(E_x,E_{\gamma })}{{\sigma }_{P_{\exp }}(E_x,E_{\gamma })}\right)}^2.$$

$$\rho (S_n)=\frac{1}{D_0}\frac{2}{g(E_x,J_{\mathrm{t}}+\text{nicefrac}12)+g(E_x,J_{\mathrm{t}}-\text{nicefrac}12)},$$

$${\rho }_{\mathrm{CT}}(E_x)=\frac{1}{T}\exp \frac{E_x-E_0}{T},$$

$$⟨{\Gamma }_{\gamma \ell }(S_n)⟩$$

$$\displaystyle\times\sum_{J_{i}}\sum_{J_{f}=|J_{i}-1|}^{J_{i}+1}g(S_{n}-E_{\gamma},J_{f})\Biggr{]},$$

$$\mathcal{L}(\mathbf{\theta })=\prod _i{\mathcal{L}}_i(\mathbf{\theta }).$$

$$\ln {\mathcal{L}}_i(\mathbf{\theta })$$

$${\chi }_{\mathrm{discrete}}^2$$

$${\chi }_{D_0}^2$$

$$+\sum _j{\left(\frac{{\rho }_{j,\mathrm{CT}}-{\rho }_{j,\mathrm{Oslo}}(\mathbf{\theta })}{{\sigma }_j}\right)}^2,$$

$${\chi }_{⟨{\Gamma }_{\gamma }⟩}^2$$

$$K_i=\sum _j\ln \frac{1}{\sqrt{2\pi {\sigma }_j}}=\sum _j\ln \frac{1}{\sqrt{2\pi \times 0.3{\rho }_{j,\mathrm{Oslo}}(\mathbf{\theta })}}.$$

$$p(E_{\gamma }|E_{\gamma }^{\prime }),$$

$$F(E_{\gamma })=\int p(E_{\gamma }|E_{\gamma }^{\prime })U(E_{\gamma }^{\prime })\mathrm{d}{E}_{\gamma }^{\prime }.$$

$$\vec{F}=P\vec{U},$$

$$AG(E_{\gamma }{)}_{E_x}$$

$$+\sum _{E_x^{\prime }<E_x}w(E_x^{\prime }{)}_{E_x}\underset{⟨AG(E_{\gamma }{)}_{E_x^{\prime }}⟩}{\underbrace{n(E_x^{\prime }{)}_{E_x}AG(E_{\gamma }{)}_{E_x^{\prime }}}}.$$

$$n(E_x^{\prime }{)}_{E_x}=\frac{S(E_x)}{S(E_x^{\prime })}.$$

$$n(E_x^{\prime }{)}_{E_x}=\frac{⟨M(E_x^{\prime })⟩N(E_x)}{⟨M(E_x)⟩N(E_x^{\prime })},$$

$$⟨M(E_x)⟩=\frac{E_x}{⟨E_{\gamma }⟩},$$

$$w(E_x^{\prime }{)}_{E_x}=\frac{FG(E_x-E_x^{\prime }{)}_{E_x}}{{\sum }_{E_{\gamma }^{\prime }}FG(E_{\gamma }^{\prime }{)}_{E_x}}.$$

$$FG(E_{\gamma }{)}_{E_x}=AG(E_{\gamma }{)}_{E_x}$$

$$-\sum _{E_x^{\prime }<E_x}n(E_x^{\prime }{)}_{E_x}\frac{FG(E_x-E_x^{\prime }{)}_{E_x}}{{\sum }_{E_{\gamma }^{\prime }}FG(E_{\gamma }^{\prime }{)}_{E_x}}AG(E_{\gamma }{)}_{E_x^{\prime }}.$$

$$\sigma (\alpha ,f)={\sigma }_{\mathrm{CN}}(\alpha )G_{\mathrm{CN}}(i,f),$$

$$G_{\mathrm{CN}}=\frac{{\Gamma }_{\gamma }(i\to f)}{{\Gamma }_{\gamma }},$$

$$n(i,f)\mathrm{d}{E}_f$$

$$={\sigma }_{\mathrm{CN}}(\alpha )\sum _{XL}\frac{⟨{\Gamma }_{\gamma }(i\to f)⟩}{{\Gamma }_{\gamma }^{(L)}}{\rho }_{\mathrm{avail}}(E_f),$$

Full text

**PROGRAM SUMMARY

*Program Title: OMpy [1]

Code Ocean Capsule: OMpy [2]

Licensing provisions: GPLv3

Programming language: Python, Cython

Reference of previous version: oslo-method-software

Does the new version supersede the previous version?: Yes

Reasons for the new version: Facilitate modular program flow and reproducible results in a transparent and well-documented code base. Updated uncertainty quantification: formerly a stage-wise normalization without built-in uncertainty propagation.

Summary of revisions: Complete reimplementation; ensemble based uncertainty quantification throughout whole method; fitting based on well-tested external libraries; corrections for the normalization procedure; auto-documentation with Sphinx

Nature of problem: Extraction of the nuclear level density and average $\gamma$-ray strength function from a set of excitation-energy tagged $\gamma$-ray spectra including the quantification of uncertainties and correlations of the results.

Solution method: The level density and $\gamma$-ray strength function can be obtained simultaneously using a set of analysis techniques called the Oslo method. To propagate the uncertainty from the counting statistics, we analyze an ensemble of perturbed spectra, which are created based on the experimental input. One obtains a set of level densities and $\gamma$-ray strength functions for each realization from a fit process. The fitting metric ($\chi^{2}$) is degenerate, but the degeneracy is removed by a simultaneous normalization of the level density and $\gamma$-ray strength function to external data, such that all correlations are preserved.

---

††journal: Computer Physics Communicationssharedfirstsharedfirstfootnotetext: Shared first authors.

1 Introduction

One long-standing challenge in nuclear physics is to precisely determine nuclear properties at excitation energies above the discrete region and up to the particle threshold(s). This region is often referred to as the quasicontinuum and represents an excitation-energy range where the quantum levels are very closely spaced. That leads to a significant degree of mixing (complexity) of their wave functions, but they are still not fully overlapping as in the continuum region. For the quasicontinuum, it is fruitful to introduce average quantities to describe the excited nucleus: instead of specific levels, the level density $\rho(E_{x})$ as a function of excitation energy $E_{x}$ is used, and instead of specific reduced transitions strengths $B(XL)$ between a given initial and final state, the average decay strength represented by the $\gamma$-ray strength function ($\gamma$SF) is applied.

In addition to their key role in describing fundamental nuclear properties, both the level density and the $\gamma$-ray strength function are vital components for calculating cross sections and reaction rates for explaining the nucleosynthesis of heavy elements in astrophysics [3, 4]. The ability to calculate cross sections may also help in the design of next generation nuclear reactors, where direct cross section measurements are missing or have too high uncertainties for the given application [5, 6, 7, 8].

The Oslo method [9, 10, 11, 12] allows for extracting the level density $\rho$ and the $\gamma$SF simultaneously from a data set of particle $\gamma$-ray coincidences. It has been implemented in a collection of programs known as the oslo-method-software [13], and the analysis has been successfully applied to a range of nuclei in widely different mass regions [14, 15, 16, 17]. However, the Oslo method consists of several non-linear steps. This makes an analytical propagation of statistical and systematic uncertainties very difficult, and thus hampers a reliable uncertainty quantification for the final results. The statistical uncertainty propagation from unfolding the $\gamma$ spectra and the determination of the primary $\gamma$-ray distribution has so far not been addressed in a fully rigorous way. In lieu of this, an approximate uncertainty estimation has been used, which is described in Ref. [12]. Moreover, uncertainties related to the absolute normalization of the level density $\rho$ and $\gamma$SF have been discussed in Ref. [18], but there was no automatized way to include this in the final results. Approximate ways to include normalization uncertainties have been suggested and used, see e.g. Ref. [19], but they do not account for correlations between parameters as they were not available within the oslo-method-software.

In this work, we approach the problem of uncertainty propagation using a Monte Carlo technique. By generating an ensemble of perturbed input spectra, distributed according to the experimental uncertainties, and propagating each ensemble member through the Oslo method, we can gauge the impact of the count statistics on the final results. The resulting level density and $\gamma$-ray strength function have to be normalized to external data. We have implemented a new simultaneous normalization for both quantities, thus preserving all correlations between them.

In the following, we present OMpy, the new implementation of the Oslo method. We discuss the various steps of the method and present our new uncertainty propagation and normalization routine. The code is tested by the analysis of a synthetic spectrum. The capability of the new method is illustrated by applying it to experimental data and a comparison to the previous implementation.

2 Usage of OMpy

OMpy is designed to facilitate a more complete uncertainty propagation through the whole Oslo method and at the same time simplify the user interface and enhance the reproducibility of the analysis. This section will focus on the latter two goals. OMpy, as well as all the Jupyter notebooks and datasets used to create the figures in this article, is available online [1, 2].

The documentation of the interface of OMpy is created automatically from its source code with the sphinx-automodapi package and is available from ompy.readthedocs.io. A detailed example of the usage is provided in the gettingstarted Jupyter notebook. Taking a step beyond just listing the version number, we have set up a docker container to ensure also the reproducibility of the software environment in which the user runs the analysis. The notebooks can be run interactively online without installation through Binder [20], although a limit on the available computation power can lead to considerable computation times for cases with large ensembles. For larger calculations we recommend the Code Ocean [DOI follows with publication] capsule.

As most operations within the Oslo method will require working on binned quantities like $\gamma$-ray spectra or level densities, or a collection of spectra into matrices, we have created the Vector and Matrix classes. These store the count data as one- or two-dimensional NumPy [21, 22] arrays, together with an array giving the energy calibration for each axis. The classes also contain a collection of convenience functions e.g. for saving and loading files to disk, rebinning and plotting.

3 The Oslo method

The starting point for the Oslo method is an $E_{x}$--$E_{\gamma}$ coincidence matrix, i.e. a set of $\gamma$-ray spectra each stemming from an identified initial excitation energy $E_{x}$. A standard way to construct this input matrix is from coincidence measurements of $\gamma$ rays and charged ejectiles following inelastic scattering or transfer reactions with light ions. An array of $\gamma$-ray detectors measures the energy $E_{\gamma}$ of the emitted $\gamma$s, while a particle telescope determines the excitation energy $E_{x}$ using the energy deposited from the outgoing charged particles. (For a detailed description, see e.g. Ref. [18] and references therein.) Recently, other methods have been developed that obtain the coincidence matrix from $\gamma$ rays in $\beta$-decay [23] or in an inverse kinematics setup [24]. An example of a coincidence matrix for 164Dy from a standard Oslo method experiment [25, 15] at the Oslo Cyclotron Laboratory is shown in panel (a) of Fig. 1.

The first step of the Oslo method is to unfold, i.e., deconvolute the $\gamma$-ray spectra for each excitation energy to compensate for the detector response (Compton scattering, pair production, etc.). This is implemented in the Unfolder class using an iterative unfolding method described in Ref. [26]. We reiterate the main points of the procedure in A. There we also describe a routine to select the best iteration, which has already been implemented in the oslo-method-software but not yet published elsewhere. The unfolded 164Dy spectrum is shown in Fig. 1 (b).

The second step is the determination of the first-generation, or primary, $\gamma$-ray spectrum for each excitation energy. In the FirstGeneration class an iterative procedure is applied as described in Ref. [27]. We recapitulate the procedure in B, including a small addition to minimize fluctuations in higher order iterations. The resulting first-generation $\gamma$-ray matrix is shown in panel (c) of Fig. 1. The main assumption of the first-generation method is that the $\gamma$-ray spectra following the populations of levels within the excitation energy bin $E_{x}$ by the inelastic or transfer reaction are the same as when this excitation energy bin is populated via $\gamma$-ray decay from higher lying levels [18]. Although this is plausible at first sight, it requires that the spin-parity distribution of the populated levels is approximately independent of the excitation energy [18, 28].

The next step of the Oslo method consists of fitting the first-generation spectra by a product of two one-dimensional functions, namely the nuclear level density $\rho(E_{x})$ and the $\gamma$-ray transmission coefficient $\mathcal{T}(E_{\gamma})$. The method relies on the relation

[TABLE]

where we look at deexcitations from an initial excitation energy bin $E_{x}$ to the final bin $E_{f}$, such that the $\gamma$-ray energy is $E_{\gamma}=E_{x}-E_{f}$. Here, $P(E_{x},E_{\gamma})$ is a matrix of first-generation spectra $FG(E_{\gamma})_{E_{x}}$ normalized to unity for each $E_{x}$ bin.333Note that we follow the standard notation for the Oslo method, where $P(E_{x},E_{\gamma})$ is the conditional probability $p(E_{\gamma}|E_{x})$ for the first $\gamma$-ray transition with energy $E_{\gamma}$ to come from an initial excitation energy $E_{x}$. Furthermore, if we assume that the $\gamma$ decay at high $E_{x}$ is dominated by dipole radiation, the transmission coefficient is related to the dipole $\gamma$-ray strength function $f(E_{\gamma})$ (or $\gamma$SF) by the relation

[TABLE]

A derivation of Eq. (1) is shown in C, where the main assumptions underpinning this decomposition are:

The compound nucleus picture: We assume that the $\gamma$ decay from an excited nuclear state is independent of how the excited state was formed. This goes back to Bohr’s theory for compound nuclei [29] and is supported by many experiments [30, 31, 32, 33, 34, 35, 15, 36, 37].

• 2.

Dominance of dipole radiation: It is assumed that the decay is dominated by dipole radiation. This is strongly supported by data and theory, see e.g. Refs. [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 37, 48].

• 3.

The generalized Brink-Axel hypothesis: The $\gamma$-ray strength function $f(E_{\gamma})$ is independent of the initial and final states, i.e., it is the same for excitations and decays between any initial and final states that are separated by the energy $E_{\gamma}$ [49, 50, 51, 52, 53, 54, 41, 55].

• 4.

The population cross-section: To obtain the ‘‘summed’’ level density $\rho(E_{f})=\sum_{J}\rho(E_{f},J)$ in Eq. (1) the population cross-section has to be approximately proportional to the intrinsic spin distribution $g(E_{x},J)$ (and we assume parity equilibration, i.e. $\rho(E_{x},J,\pi=+)\approx\rho(E_{x},J,\pi=-)$). When the employed reaction is very spin selective, like in $\beta$-decay, a weighted sum of the level densities is obtained instead, see Eq. (40).

Finally, we have to normalize the nuclear level density $\rho(E_{x})$ and $\gamma$-ray strength function $f(E_{\gamma})$ (or equivalently, the $\gamma$-ray transmission coefficient $\mathcal{T}(E_{\gamma})$). It was shown by Schiller et al. [12] that the solution to Eq. (1) is invariant under a Lie group $G$ of transformations by three continuous real valued parameters $A,B$ and $\alpha$:

[TABLE]

However, we stress that the degrees of freedom are limited to those given by $G$ --- i.e., all shape features of $\rho$ and $f$ beyond the exponential prefactor given by $G$ are uniquely determined by the fit. It is important to note that the $\alpha$ parameter, which influences the slope (in a log plot) of the functions, is common to $\rho$ and $f$. Hence, their normalizations are coupled together in the Oslo method. To obtain the physical solution, the level density and $\gamma$SF need to be normalized to auxiliary data. Typically, one uses $s$-wave resonance spacings, $D_{0}$, from neutron capture experiments as well as discrete levels to fix the level density normalization, and augment this by average total radiative width data, $\langle\Gamma_{\gamma}\rangle$, to normalize the $\gamma$SF. This will be discussed in more detail in the next sections.

4 Uncertainty propagation by ensemble

We use an approach based on the Monte Carlo (MC) technique to estimate the statistical uncertainties in the Oslo method by creating an ensemble of randomly perturbed copies of the data set under study. To illustrate the method, we have chosen an experimental data set for 164Dy. The data was obtained in Ref. [25] and we will compare our results to the reanalysis published by Renstrøm et al. [15].

The random variables are the experimental number of counts in each bin $i$ in the raw $E_{x}$--$E_{\gamma}$ coincidence matrix $R$. We assume that they are independent and follow a Poisson distribution with parameter $\lambda_{i}$. The Poisson distribution $\mathcal{P}_{\lambda}$ is given as

[TABLE]

We take the number of counts $k_{i}$ in bin $i$ of $R$ as an estimate for the Poisson parameter $\lambda_{i}$. Note that it is an unbiased estimator for $\lambda_{i}$, since the expectation value $\langle k\rangle=\lambda$. To generate a member matrix $R_{l}$ of the MC ensemble, the counts in each bin $i$ are replaced by a random draw from the distribution $\mathcal{P}_{k_{i}}$. By this procedure, we obtain $N_{\mathrm{ens}}$ matrices representing different realizations of the experiment. Defining $\vec{r}_{i}=(r_{i}^{(1)},r_{i}^{(2)},\dots,r_{i}^{(N_{\mathrm{ens}})})^{T}$ as the vector of all $N_{\mathrm{ens}}$ realizations $m$ of bin $i$ we can calculate the sample mean $\langle\vec{r}_{i}\rangle$,

[TABLE]

and standard deviation $\sigma_{i}$,

[TABLE]

Of course, in the case of the raw matrix $R$, the standard deviation is trivial because it is given by the Poisson distribution ($\sigma=\sqrt{\lambda}$). But the technique also allows us to estimate the standard deviation at later stages in the Oslo method --- after unfolding, after the first-generation method and even after fitting the level density and $\gamma$-ray transmission coefficient. In Fig. 1 we show the relative standard deviations $\sigma_{\mathrm{rel},i}=\sigma_{i}/r_{i}$ in the raw (d), unfolded (e) and first-generation (f) matrices of the 164Dy dataset based on $N_{\mathrm{ens}}=50$ ensemble members.

It should be noted that we usually analyze spectra with a (possibly time-dependent) background. In this case we measure two raw spectra, the total and the background spectra, which independently follow a Poisson distribution. In OMpy, both spectra can be read in and perturbed independently. The background subtracted spectra $R^{\prime}$ are then generated for each realization. When the number of background counts are large enough, this may lead to bins in $R^{\prime}$ with negative number of counts. With the current default, these are removed before further processing at the cost of a potential bias towards higher level densities $\rho$ and strength functions $f$ (see the discussion in Sec. 7). The generation of the all ensembles matrices (including application of the unfolding and first-generation method) is handled by the Ensemble class.

5 Decomposition into level density and $\gamma$-ray

strength function

With the first-generation matrices and their standard deviation at hand, we may proceed with the fitting of $\rho$ and $\mathcal{T}$ for each ensemble member. First, we select a suitable bin size $\Delta E$, typically 100--300 keV depending on the detector resolution, and rebin the first-generation matrix along both the $E_{x}$ and $E_{\gamma}$ axes to this bin size. The matrix of experimental decay probabilities $P_{\mathrm{exp}}(E_{x},E_{\gamma})$ is obtained by normalizing the spectrum in each $E_{x}$ bin to unity. For the fit of $\rho$ and $\mathcal{T}$, we take the function value in each bin as a free parameter. For a given pair of trial functions $(\rho,\mathcal{T})$, the corresponding matrix $P_{\mathrm{fit}}(E_{x},E_{\gamma})$ is constructed by

[TABLE]

where $C_{E_{x}}$ is a normalization coefficient so that $\sum_{E_{\gamma}}P_{\mathrm{fit}}(E_{x},E_{\gamma})=1\,\,\forall E_{x}$. We fit $P_{\mathrm{fit}}$ by a $\chi^{2}$ minimization approach, minimizing the weighted sum-of-squared errors

[TABLE]

It is important to use a weighted sum rather than simply a sum of the residuals, to suppress the influence of bins with large uncertainties. This in turn makes uncertainty estimation important. As already mentioned, a shortcoming of the original Oslo method implementation [12] in the oslo-method-software has been the estimation of the uncertainty $\sigma_{P_{\mathrm{exp}}}(E_{x},E_{\gamma})$ in the denominator of the $\chi^{2}$ fit. Due to the lack of a complete statistical uncertainty propagation, one has had to resort to an approximate uncertainty estimation. It was based on a Monte Carlo scheme similar to the present work, but where only the first-generation spectrum is perturbed as if each entry was direct count data. This is discussed in detail in Ref. [12]. In OMpy, we have access to a proper uncertainty matrix $\sigma_{P_{\mathrm{exp}}}$. We checked that most bins of the first-generation matrices approximately follow a normal distribution. However, they are distributed with a larger and varying standard deviation as compared to what one would have received if the first generation entries were count data directly and followed the expectation value $\langle k\rangle=\lambda$. The approximate adherence to the normal distribution justifies the use of a $\chi^{2}$ minimization as a likelihood maximization.

The $\chi^{2}$ minimization is carried out by numerical minimization in the Extractor class. This is a different implementation than in the oslo-method-software, where the minimum is found by iteratively solving a set of equations to obtain a solution satisfying $\partial\chi^{2}/\partial\rho=0$, $\partial\chi^{2}/\partial\mathcal{T}=0$ for each bin of $\rho$ and $\mathcal{T}$ [12]. After testing several off-the-shelf minimizers, we have found that the modified Powell’s method in the SciPy package works well [56, 57].

The normalized first-generation matrix $P(E_{x},E_{\gamma})$ is compared to the fitted matrix in Fig. 2 for one ensemble member of the 164Dy dataset. In Fig. 3 we show the corresponding level density $\rho$ and $\gamma$-ray strength function $f$, where the latter is obtained from $\mathcal{T}$ using Eq. (2).

5.1 Testing the sensitivity on the initial values

The presently used minimization routine requires an initial guess for the trial functions $(\rho,\mathcal{T})$. In principle, the choice of the trial functions may have a significant effect on the results if the minimizer is prone to get caught in local minima. Ref. [12] proposes to set the initial $\rho$ to a constant, arbitrarily chosen as 1, for all excitation energies $E_{x}$, and $\mathcal{T}(E_{\gamma})$ as the projection of $P_{\mathrm{exp}}(E_{x},E_{\gamma})$ on the $E_{\gamma}$-axis. In OMpy, we have implemented several choices to test the stability of the solutions towards other initial guesses. For the cases shown here, no significant impact on the final results was observed.

Our default choice for the initial guess on $\rho$ has been motivated by a long lasting discussion on whether the level density $\rho$ follows the back-shifted Fermi gas (BSFG) equation [58, 59] or a constant temperature (CT) model [60, 59] below the neutron separation energy $S_{n}$ , see Ref. [61, 62] and references therein. In recent years, the results of the Oslo method have strongly suggested a close-to CT-behavior. As this is equivalent to the initial guess for $\rho$ in Ref. [12] (after a transformation $G$ (see Eq. (3)), we set the default initial guess to a BSFG-like solution. More specifically, we draw from a uniform distribution centered around a BSFG-like initial guess for each ensemble member. Fig. 3 demonstrates that the resulting fit is still best described by the CT-model, but it is now very unlikely that this is due to a failure of the minimizer. The initial guess for $\mathcal{T}$ is still chosen as in Ref. [12], but with the same randomization with a uniform distribution as for the level density $\rho$.

Besides the two alternatives named above, we have also tested a rather exotic initial guess for $\rho$ given by a quadratic function with a negative coefficient for the second-degree. This contrasts any expectation that the level density $\rho$ increases as a function of the excitation energy $E_{x}$. However, even for this choice, the solution is stable.

5.2 Uncertainty estimation

Given the degeneracy of the $\chi^{2}$, it is nontrivial to estimate the uncertainty in the solutions of $\rho$ and $\mathcal{T}$. The oslo-method-software implements the approach of Ref. [12], where the uncertainty is estimated from the standard deviation of the solutions ($\rho^{(m)}$, $\mathcal{T}^{(m)}$) for each realization of the first-generation matrix $P_{\mathrm{exp}}^{(m)}$. This may lead to erroneously high uncertainties, as any transformation $G$ gives an equally good solution to a given $P_{\mathrm{exp}}^{(m)}$. To illustrate this, we could imagine that the $P_{\mathrm{exp}}^{(m)}$’s only differ by a very small noise term. Due to the noise term the solutions will not be identical. In addition, we recapitulate that any given $P_{\mathrm{exp}}^{(m)}$ can be equally well fit by $\rho$ as with any allowed transformation $G$ such as 10$\rho$ or 100$\rho$. Instead of a negligible standard deviation one receives a standard deviation solely based on the degeneracy of the solutions. So far, we have observed that the minimization results are more stable as one might expect for the scenario outlined here (provided that the initial guess is not randomized). This is probably due to the way the minimizer explores the parameter space. Nevertheless, for the standard usage of OMpy we provide the functionality to estimate the uncertainty of solutions $\rho$ and $\mathcal{T}$ after normalization of each ensemble of ($\rho^{(m)}$, $\mathcal{T}^{(m)}$). This will be explained in more detail in the next section.

5.3 Fit range

Using the Oslo method we have to restrict the fit range for the $P(E_{x},E_{\gamma})$ matrix. The area below $E_{x}^{\mathrm{min}}$ exhibits discrete transitions, thus does not adhere to the statistical nature of the $\gamma$-ray strength function $f$ and is therefore excluded. To remain selective on the $\gamma$-decay channel, we can only use $P(E_{x},E_{\gamma})$ up to a maximum excitation energy $E_{x}^{\mathrm{max}}$ around the neutron separation energy $S_{n}$. Finally, we also constrain the minimum $\gamma$-ray energy to $E_{\gamma}^{\mathrm{min}}$, which usually attributed to a deficiency of the unfolding or first-generation method for low $\gamma$-ray energies. The limits are highlighted in Fig. 1 (c) and only this valid region of $P(E_{x},E_{\gamma})$ is sent to the minimizer.

6 Normalization of $\rho$ and $\gamma$SF

At first glance, the results bear little resemblance to a level density or $\gamma$SF. That is because the fit has not yet been normalized. Thus, the solution shown in Fig. 3 is just one of an infinite set of solutions to the fit. In this section we will discuss how auxiliary data can be used to find the transformation $G$ that gives the physical solutions.

6.1 Auxiliary experimental data

For the level density $\rho$, there are often two different types of auxiliary datasets available. At low energies, the discrete levels are known from spectroscopy. They can be compared to the fitted level density from the Oslo method after applying the same binning. One can also account for the detector resolution by applying e.g.  a Gaussian smoothing to the histogram. At higher excitation energies, typically $\sim$1-3 MeV, the spectroscopy data fails to resolve all levels [61]. The user will thus have to set a sensible region in the low energy regime for the comparison.

The second piece of information stems from neutron resonance experiments, e.g. the average $s$-wave resonance spacings $D_{0}$. They provide information about the level density $\rho(E_{x}\,{=}\,S_{\mathrm{n}},J_{\mathrm{t}}\pm\nicefrac{{1}}{{2}},\pi_{\mathrm{t}})$ at the neutron separation energy $S_{n}$, where $J_{\mathrm{t}}$ and $\pi_{\mathrm{t}}$ are the ground-state spin and parity of the target nucleus, i.e.  the $A-1$ nucleus [61]. With the Oslo method we obtain the total level density $\rho(E_{x})=\sum_{J,\pi}\rho(E_{x},J,\pi)$. If one knows the fraction of $J_{\mathrm{t}}\pm\nicefrac{{1}}{{2}},\pi_{\mathrm{t}}$ levels, one can estimate $\rho(S_{n})$ by from $D_{0}$ by [12]

[TABLE]

where $g(E_{x},J,\pi)$ is the spin-parity distribution of the nucleus at $E_{x}$ and the factor of 2 comes from the assumption of equiparity, i.e.  $g(E_{x},J,\pi)=g(E_{x},J)/2$. Note that the $J_{\mathrm{t}}-\nicefrac{{1}}{{2}}$ term vanishes for $J_{\mathrm{t}}=0$. For the spin-parity distribution $g(E_{\mathrm{x}},I)$ it is common to use a form proposed in Ref. [63, Eq. (3.29)]; however the exact parametrization is a major source of systematic uncertainties in the normalization, thus several suggestions are implemented in OMpy.

The usage of $\rho(S_{n})$ for the normalization is further complicated by restricted fit-regions for $P(E_{x},E_{\gamma})$ (see Sec. 5.3). These limit the extraction of the level density $\rho$ up to $E_{x}^{\mathrm{max}}$--$E_{\gamma}^{\mathrm{min}}$, which is often about 1-3 MeV below $S_{n}$. Consequently, we cannot directly normalize the fitted level density $\rho$ to $\rho(Sn)$ obtained from $D_{0}$. To utilize this information, we have to extrapolate $\rho$ and compare the extrapolation at $S_{n}$. The exact form of the extrapolation is another systematic uncertainty. Generally, the constant temperature (CT) model [63] fits well with the level density data obtained from the Oslo method [64],

[TABLE]

where the temperature $T$ and the energy shift $E_{0}$ are free parameters of the model. We obtain them by a fit to $\rho$ in a suitable energy range.

The convolution of the level density $\rho$ and $\gamma$-ray strength function $f$ (or equivalently $\mathcal{T}$) can be further constrained by the average total radiative width $\langle\Gamma_{\gamma}(S_{n})\rangle$ from neutron-capture experiments (see e.g. Ref. [61]) on a target nucleus with ground-state spin $J_{t}$ using the following equation:

[TABLE]

where the first sum runs over all possible residual nucleus spins $J_{i}$, i.e.  from $\min|J_{t}\pm\nicefrac{{1}}{{2}}\pm\ell|$ to $J_{t}+\nicefrac{{1}}{{2}}+\ell$ and the second sum runs over all final spins $J_{f}$ accessible with dipole radiation starting from a given $J_{i}$. Often only s-wave information is available, corresponding to $\ell=0$ in the notation above, and the measurements are performed with low energy neutrons, such that $E_{x}\approx S_{n}$; for brevity it is then common to write $\langle\Gamma_{\gamma}\rangle$ for $\langle\Gamma_{\gamma\ell}(S_{n})\rangle$. A derivation of Eq. 6.1 is given in D. Currently, we approximate $\langle\Gamma_{\gamma}\rangle$ using the integral involving the level density $\rho$ for all excitation energies $E_{x}$, although the integral at low energies can be replaced by a sum over decays to the known discrete levels for more precise calculations. Note that we also have to extrapolate the $\gamma$-ray strength function to be able to use this equation. Due to its shape, a log-linear function is often fitted to the results.

6.2 Likelihood

In the following, we will define the likelihood $\mathcal{L}(\bm{\theta})$ that is used to find the proper normalizations. Let $\mathcal{L}_{i}(\bm{\theta})$ denote the likelihood for a given solution ($\rho$, $\mathcal{T}$) of Eq. (7) to match the normalization information $i$ (i.e.  discrete levels, $D_{0}$, $\langle\Gamma_{\gamma}\rangle$) after the transformation $G$ with the parameters $\bm{\theta}=(A,B,\alpha)$. Due to the extrapolation of the level density $\rho$ mentioned above, we have to introduce two nuisance parameters $T$ and $E_{0}$, so we extend $\bm{\theta}$ to include $(A,B,\alpha,T,E_{0})$. To reduce the computational complexity, we extrapolate the $\gamma$-ray strength function by the best-fit values for a given set of transformations $(\alpha,B)$. The total likelihood $\mathcal{L}(\bm{\theta})$ is then given by

[TABLE]

We assume that the experimental $D_{0}$ and $\langle\Gamma_{\gamma}\rangle$ data are normal distributed, thus maximizing the log-likelihood is equivalent to minimizing a sum of $\chi^{2}_{i}$’s. Note that the measurement of the discrete levels is of course not stochastic; however, the count data we use to determine $\rho$ for the comparison is. More explicitly, we have

[TABLE]

where $K_{i}=\ln(1/(2\pi\sigma_{i}))$ is a constant as long as the standard deviation(s) $\sigma_{i}$ does not depend on $\bm{\theta}$. The subscript exp denotes the experimental data. The sums in Eq. (14) and (6.2) run over all data points used in the evaluation, and for $\chi^{2}_{D_{0}}$ we invert Eq. (9) to obtain $D_{0}$ from the level density (extrapolated with the CT model). The second term of Eq. (6.2) arises due to the fit of the nuisance parameters of the CT model ($T$, $E_{0}$) to the Oslo method data, here labeled $\rho_{\mathrm{Oslo}}$.

As discussed in Sec. 5.2, the degeneracy of the solutions $\rho$ and $f$ prevents us from directly inferring their parameter uncertainties in the fit of the first generation matrix $P_{\mathrm{exp}}$. However, clearly the data points of $\rho$ and $\mathcal{T}$ have a statistical uncertainty that should propagate to inform the posterior distribution of $\bm{\theta}$. We choose to model this by setting $\sigma_{j}$ somewhat arbitrarily to a relative uncertainty of 30%. We propose to test the implications in a future work e.g. by comparing inferred $D_{0}$’s for datasets where $D_{0}$ is known but on purpose not included in the $\chi^{2}$ fit of Eq. (6.2). Note that with this specification of the standard deviations $\sigma_{j}$, the $K_{i}$ become $\bm{\theta}$-dependent:

[TABLE]

The likelihood can easily be extended if other information shall be taken into account. In several recent works on the Oslo method, experimental data on $\langle\Gamma_{\gamma}\rangle$ was not available, but roughly estimated from systematics (see e.g. Ref. [65, 19]). With the new possibilities of OMpy we would instead recommend to constrain $\alpha$ and $B$ by adding a term to the total likelihood that describes the match with other measured strength function data (usually above $S_{n}$).

6.3 Implementation

We sample this likelihood with the Bayesian nested sampling algorithm MultiNest [66, 67] using the PyMultiNest module [68]. For a more efficient calculation, we first find an approximate solution (more accurately the maximum-likelihood estimator) $\bm{\hat{\theta}}$ with the differential evolution minimizer of SciPy [56, 69]. This is by default used to create weakly informative priors for $A$, $B$ and $\alpha$ and $T$. For $A$ and $B$ we use a normal distribution truncated at 0 (negative values of $\rho$ or $f$ are not meaningful), a default mean $\mu$ given by $\bm{\hat{\theta}}$ and a broad width of $10\mu$. For $\alpha$ and the nuisance parameter $T$ (entering through the level density extrapolation model) we use log-uniform priors spanning one order of magnitude around $\bm{\hat{\theta}}$. For the second nuisance parameter $E_{0}$ we choose a normal distribution with mean 0 and width 5 MeV, which is truncated below (above) -5 and 5 MeV, respectively. The latter choice is well justified regarding the range of reported values of $E_{0}$ in Ref. [70].

This simultaneous normalization is implemented in the NormalizerSimultan class, which relies on composition of the NormalizerNLD and NormalizerGSF classes handling the normalization of the level density $\rho$ and $\gamma$-ray strength function $f$. To facilitate a comparison with the oslo-method-software calculations, one can also run a sequential normalization first using Eq. (14) and (6.2) through NormalizerNLD and then the resulting $\rho$ as input to the normalization through NormalizerGSF, see Eq. (6.1). It should be stressed though that the normalization through the likelihood calculations in OMpy will still differ from the approach taken in the oslo-method-software. The latter allowed only to receive best-fit estimates of the transformation parameters $A$, $B$ and $\alpha$. Any subsequent uncertainty estimation due to the normalization itself was up to the users. We also note that an advantage of the simultaneous approach is that one obtains the correlations between $A$, $B$ and $\alpha$ such that uncertainties in the normalization of $\rho$ directly propagate to the estimation transformation parameters for $f$.

From the MultiNest fit we obtain the posterior probability distribution for the parameters $\bm{\theta}$, given as (equally-weighted) samples $\bm{\theta}_{i}$.444For brevity, we drop the index $m$ on $\bm{\theta}$, but the normalization is performed for each realization $m$ individually. The normalization uncertainty for the solution ($\rho^{(m)}$ , $f^{(m)}$) of the realization $m$ can then be mapped out by creating a normalized sample $(\rho^{(m)},f^{(m)})_{i}$ for each $\bm{\theta}_{i}$ using Eq. (3). By repeating this procedure for all realizations $m$ of the ensemble, we also recover the uncertainty due to the counting statistics. This is performed in the EnsembleNormalizer class.

7 Systematic uncertainties

The previous sections described the necessary tools to evaluate the statistical uncertainties due to the counting statistics and the normalization procedure. It is important to keep in mind that there are also systematic uncertainties linked to the analysis, which are summarized below:

1. a)

Removal of negative counts in raw matrix: When subtracting the background from the raw matrix, we often receive a matrix with negative counts in some bins. This is clearly linked to the Poisson statistics in regions with a low signal to background ratio. If one simply removes the negative counts, one potentially biases the mean of the level density and $\gamma$-ray strength points derived from these bins. It was previously observed that negative counts in the raw matrix can cause technical challenges in the currently implemented unfolding method, with some bins obtaining extreme negative values after several iterations. For the background subtraction, it might be a more reasonable fix to redistribute the negative counts to bins within the resolution. This is implemented by Matrix.fill_negative as an alternative to the default method, Matrix.remove_negative, which removes all negative counts. For the 164Dy, there is a high signal to background ratio, thus the background subtraction does not have this problem. In cases with a worse background ratio, both methods could be compared. If they result in significant differences, further analysis is needed to find the optimal procedure. 2. b)

Removal of negative counts in unfolded and first-generation matrix: The unfolding and first-generation methods can result in negative counts which can not be linked to the Poisson statistics any longer. It is thus not clear whether it is better to keep the negative counts or to redistribute them in the fashion described above. Again, the problem was negligible in the 164Dy dataset, but has to be treated carefully if the methods lead to more bins with such a behavior. To retrieve the matrices before removal of the negative counts, the user can simply replace or deactivate the remove_negative methods of the Unfolder and FirstGeneration classes. 3. c)

Unfolding method: There are two main sources of systematic uncertainties, the iterative unfolding method itself, and the detector response functions. The latter can be gauged by obtaining an ensemble of different detector response functions that capture the breadth of physically reasonable configurations (e.g. auxiliary software such as GEANT4 [71]). Each member of the raw matrices $R^{(m)}$ is then unfolded with one (or each) of the different detector response functions. The former uncertainty is more difficult to quantify. In this special case, alternative methods exist already and one approach could be to implement an alternative unfolding algorithm (e.g. Ref. [72]). The oslo-method-software attempts to quantify the systematic uncertainty from unfolding and the first-generation matrix following an ad hoc numerical procedure described in Ref. [12]. A proper treatment that pays tribute to the full complexity of the problem is outside the scope of the present work. 4. d)

Population cross-section: As mentioned in Sec. 3, the first-generation method assumes that the spin-parity distribution of the populated levels $g_{\mathrm{pop}}$ is similar for the whole excitation energy range studied. At best, $g_{\mathrm{pop}}$ approximates the spin-parity distribution $g$ of the levels in the nucleus itself. This is often believed to be the case for low and mid-mass nuclei, where the beam energy is chosen such that the compound cross-section dominates over the direct cross-section. More details can be found in Refs. [18, 73] for challenges if the intrinsic spin-parity distribution $g(E_{x},J,\pi)$ is very different from the (normalized) population cross-section $g_{\mathrm{pop}}(E_{x},J,\pi)$. 5. e)

Decomposition and the Brink-Axel hypothesis: The decomposition of the first-generation matrix into the level density $\rho(E_{x})$ and $\gamma$-ray strength function $f(E_{\gamma})$ relies on a generalization of the Brink-Axel hypothesis [49, 50], where the strength function is assumed to be approximately independent of $E_{x}$, $J$ and $\pi$ (see Eq. (C)). The validity of the assumption has been tested within the Oslo method by comparison of the $\gamma$-ray strength function $f$ extracted from different initial and final excitation energy bins, see Refs. [51, 52] and references therein for further works. 6. f)

Intrinsic spin-parity distribution: Both the normalization of the level density $\rho$ at $S_{n}$ via $D_{0}$ and the absolute scaling of the $\gamma$-ray strength function via $\langle\Gamma_{\gamma}\rangle$ rely on the knowledge of the intrinsic spin-parity distribution $g(E_{x},J,\pi)$ of the nucleus, see Eq. (9) and Eq. (6.1). It is difficult to measure the spin-parity distribution $g(E_{x},J,\pi)$ in the (quasi)continuum and there are various different empirical parametrization and theoretical calculations, see e.g. Refs. [70, 74, 59, 75] and Refs. [76, 77, 78, 79] respectively. 7. g)

Lack of $D_{0}$ and $\langle\Gamma_{\gamma}\rangle$: In several recent cases where the Oslo method has been applied, experimental values of $D_{0}$ and/or $\langle\Gamma_{\gamma}\rangle$ were not available, see e.g. Refs. [65, 19, 15, 37]. In these works, $D_{0}$ and/or $\langle\Gamma_{\gamma}\rangle$ were estimated from the values of nearly nuclei. With OMpy it is now easier to adopt the normalization procedure via the likelihood $\mathcal{L}(\bm{\theta})$ in Eq. (12) instead. One can either remove terms where the experimental information is missing, allowing for a larger degeneracy of the solutions, or add different constrains, like a measure of how well the $\gamma$-ray strength function $f(E_{\gamma}$) matches other measured $\gamma$-ray strength function data, which often exists above $S_{n}$. 8. h)

Bin sizes and fit ranges: Several decisions have to be made on the bin sizes and fit ranges, for example in the normalization of the level density $\rho$ to the known level scheme, see Sec. 6.1. The user can test different sensible regions for this comparison.

In this work, only the impact of the negative counts, items a and b have been studied. In many cases it is challenging to evaluate the impact of the assumptions on the data. Whenever the uncertainty is linked to models, algorithms, parameter sets or alike, and alternative choices exists (items c, d, f, h), it is possible to utilize OMpy’s error propagation functionality. This was illustrated for the unfolding of the data in item c. More information on items e and f can be found in Ref. [18].

8 Discussion and comparison

In Fig. 4 we show the level density $\rho$ and $\gamma$-ray strength function $f$ of 164Dy resulting from the simultaneous normalization of $N_{\mathrm{ens}}=50$ realizations. Each realization $m$ is transformed with $N_{\mathrm{samp}}=100$ samples from the normalization parameters $\bm{\theta}_{i}$. The combined uncertainty of the normalization and counting statistics is visualized through the 16th, 50th, and 84th percentiles which together form the median and a 68% credible interval. Additionally, we show one randomly selected sample $(\rho^{(m)},f^{(m)})_{i}$ including its extrapolation. The results are compared to the analysis of Renstrøm et al. [15] which used the oslo-method-software, and we display both the uncertainty that is quoted due to the counting statistics and the total uncertainty, that includes the normalization. Note that we discussed in Sec. 5.2 that this split into the uncertainty due to counting statistics and due to normalization is, strictly speaking, not possible -- hence, the quoted decomposition is an approximation.

It is gratifying to see that overall, both analyses provide similar results. In Fig. 4a, the level density below $\sim 2$ MeV exhibits the same structure of bumps attributable to the discrete level structure, and at higher $E_{x}$ the curves are practically identical. Similarly, for the $\gamma$SF in Fig. 4b, the data is mostly compatible within the error-bars. The median of the results from OMpy has slightly steeper slope and the uncertainties are somewhat more evenly distributed across the whole energy range.

Some deviation between the results is expected due to the different fitting method and software implementations. Moreover, different fit-regions may lead to different estimations of the normalization parameters $\theta$, which in turn give different slopes (and absolute values) for $\rho$ and the $\gamma$SF. Usually, there will be several sensible fit-regions. With OMpy, the user could create a wrapper that loops through the different fit-regions. The mean and spread of the results can be analyzed using the same ensemble based approach as above.

Another way to verify the results of OMpy is to simulate decay data for a given $\rho$ and the $\gamma$SF and compare these to the analysis with OMpy from the simulated data. We have used the Monte-Carlo nuclear decay code RAINIER v1.5 [81, 82] to create decay data from a 164Dy-like nucleus. To create an artificial level scheme, we used the first 20 known discrete levels [61] and for higher energies created levels following the CT model and spin-distribution described in Renstrøm et al. [15]. The $\gamma$SF has been modeled with the parameters from the same publication. Then, we simulate the experiment by populating $2\times 10^{6}$ levels below $S_{n}$ and recording the decay $\gamma$s for each event. For simplicity, we have assumed here that we populate the levels proportionally to the intrinsic spin distribution of the nucleus. Finally, we use the response functions of the 164Dy experiment to convert the recorded $\gamma$ rays to a matrix of synthetically generated events that substitute the experimentally determined raw matrix in the further analysis with OMpy. The full setting file can be found with the supplementary material online.

In Figure 5 we compare the results of OMpy to the input level density $\rho$ and $\gamma$SF from RAINIER. The fitted level density $\rho$ slightly over-pronounces the structures of the discrete levels at low excitation energies. From a practical point of view this is not a significant problem, as the Oslo method is used for an analysis of the quasicontinuum region, so at energies $\gtrsim 2$ MeV. There, we observe a perfect match between the input and the fitted level density. The resulting $\gamma$SF is in very good agreement with the input $\gamma$SF. There are small deviations at the lowest and highest energies that can easily be understood. The apparent discrepancy of the input $\gamma$SF below $\sim 2$ MeV can be explained by a failure of the first-generation method due to strong populations of discrete states and is directly visible in the comparison of the unfolded to first-generation matrix (see supplementary material online). The region could have been excluded from the comparison, but we aimed for the same extraction region as for the experimental 164Dy data set. The mismatch propagates to a wrong extrapolation of the $\gamma$SF towards lower energies. The two resonance-like structures at $\sim$6 and 7 MeV over-compensate for the mismatch of the level density at the lowest excitation energies.

One of the main motivations for OMpy was to improve the uncertainty analysis in the Oslo method. Throughout the article we have highlighted several improvements to the theoretical framework for propagating the uncertainties compared to the oslo-method-software. The most fundamental source of uncertainties in the analysis is the count statistics. For each experiment a balance between the run-time costs and the possibility to gather more data has to be found. Therefore, we find it instructive to study the impact of a reduced number of counts for the 164Dy experiment. We create a new raw matrix, which we draw from a Poisson distribution with a mean of 1/10th of the original data (an arbitrary choice corresponding to 1/10th of the run actual time) and otherwise perform the same analysis as above. Both analyses are compared in Fig. 6, and overall we find a good agreement of the median values.

There is one major difference in the analysis that leads to different uncertainty estimates. In the case with the reduced number of counts we have to use a lower maximum excitation energy $E_{x}^{\mathrm{max}}$ in the fit to $P(E_{x},E_{\gamma})$. As a consequence, the extracted level density ranges only up to about 4.5 MeV instead of 6.5 MeV, increasing the uncertainty in the determination of the normalization parameters. The increased uncertainty of the slope parameter $\alpha$ is directly visible for the level density $\rho$. For the $\gamma$SF we additionally require a suitable fit to the experimental $\langle\Gamma_{\gamma}\rangle$, so a different slope $\alpha$ will also affect the absolute scaling $B$. This results in the normalized $\gamma$SFs tilting around a $\gamma$-ray energy of about 5 MeV, where the exact location of the tilting point depends on the nucleus. The normalization uncertainty thus leads to increasing relative uncertainties of the $\gamma$SF for lower $\gamma$-ray energies $E_{\gamma}$. This contrasts to the naive expectation that the relative uncertainties should increase with $E_{\gamma}$, as the number of counts that determine each $\gamma$SF bin decrease. Finally, we foresee that the updated framework for the uncertainty analysis may have a significant impact when further processing of the results from the Oslo method. An example of this is given in Renstrøm et al. [15], where the authors fit of the strength of the peak at $\sim 3$ MeV in the $\gamma$SF which assumed to be a $M1$ scissors mode.

9 Extension of OMpy

OMpy is written with modularity in mind. We want it to be as easy as possible for the user to add custom functionality and interface OMpy with other libraries. For example, in Sec. 7 we discussed that it may be of interest to try other unfolding algorithms than the one presently implemented. To achieve this, one just has to write a wrapper function with the same return structure as the callable Unfolder class. Then one provides the new wrapper to Ensemble instead of the Unfolder class and all matrices will be unfolded with the new algorithm.

It is our hope and goal that OMpy will be used as much as possible. Feedback and suggestions are very welcome. We encourage users who implement new features to share them by opening a pull request in the Github repository.

10 Conclusions and outlook

We have presented OMpy, a complete reimplementation of the Oslo method in Python. The capabilities of the code have been demonstrated by comparison with a synthetic data set modeled with the decay code RAINIER. We have also compared OMpy to the analysis for 164Dy with the previous implementation of the Oslo method, the oslo-method-software and in general find a good agreement. However, we have refined the uncertainty quantification of the analysis using an ensemble-based approach. We are now able to simultaneously take into account statistical uncertainties from the counting statistics and the normalization to external data and preserve the full correlation between the resulting level density $\rho$ and $\gamma$SF.

Acknowledgements

We would like to thank V. W. Ingeberg for a beta test and his feedback on the software. This work was supported by the Research Council of Norway under project Grants No. 263030, 262952 and 23054. A.C.L. gratefully acknowledges funding from the European Research Council, ERC-STG-2014 Grant Agreement No. 637686.

Appendix A Unfolding

Here we explain the unfolding technique presented in Ref. [26], which is used both in the original Oslo method implementation in the oslo-method-software and in OMpy. Let the detector response be modeled as a conditional probability density function

[TABLE]

encoding the probability that a $\gamma$ ray with true energy $E_{\gamma}^{\prime}$ is detected with energy $E_{\gamma}$. Given a true $\gamma$-ray spectrum $U(E_{\gamma})$, the folded spectrum $F(E_{\gamma})$, i.e. the spectrum seen by the detector, is then given by

[TABLE]

By discretizing into energy bins of width $\Delta E_{\gamma}$, it becomes a matrix equation

[TABLE]

where $P$ is the response matrix of discrete probabilities $P_{kl}=p(E_{\gamma,k}|E_{\gamma,l}^{\prime})\Delta E_{\gamma}$. The unfolding procedure amounts to solving this equation for $\vec{U}$. However, a straightforward matrix inversion is ill-advised, as it will often lead to singularities or produce large, artificial fluctuations in $\vec{U}$ [83, 84]. Instead, the approach taken in the Oslo method is to use an iterative technique that successively approximates $\vec{U}$. Letting $\vec{R}$ denote the measured spectrum, the algorithm is

Start with a trial function $\vec{U}_{0}=\vec{R}$ at iteration $i=0$ 2. 2.

Calculate the folded spectrum $\vec{F}_{i}=P\vec{U}_{i}$ 3. 3.

Update the trial function to $\vec{U}_{i+1}=\vec{U}_{i}+(\vec{R}-\vec{F}_{i})$ 4. 4.

Iterate from 2 until $\vec{F}_{i}\approx\vec{R}$.

Note that the oslo-method-software uses a custom tailored combination of the additive updating procedure of step 3, and a ratio approach, $\vec{U}_{i+1}=\vec{U}_{i}\circ(\vec{R}\oslash\vec{F}_{i})$, where the Hadamard products stand for an element wise product. We obtain equally good results by adopting only the additive updating in OMpy. The updating procedure may sometimes lead to a negative number of counts in the unfolded spectra. For negative counts close to zero, it is not clear whether these should be kept or discarded, as it is not clear whether they originate from the statistical nature of the data. In some cases, one observes some bins with large, negative counts which hints at a failure of the method. These cases should be analyzed carefully before the results are processed further. The current default is to remove the negative counts at the end of the unfolding.

A too large number of iterations does not improve the results significantly, but introduces strong fluctuations in the unfolded spectrum. After the publication of Ref. [26] a criterion for step 4 has been added to the oslo-method-software, which is also used in OMpy. A predefined number (usually around 30-200) iterations is run and the best iteration is selected based on a weighted sum over each vector element of the root-mean-square error of $\vec{F}_{i}-\vec{R}$ and the relative level of fluctuations in $\vec{U}_{i}$ compared to the fluctuations of the raw spectrum $\vec{R}$. The relative fluctuations are estimated as $\lvert(\vec{U}_{i,l}-\langle\vec{U_{i}}\rangle)/\langle\vec{U_{i}}\rangle\rvert_{1}$, where $\langle\vec{U_{i}}\rangle$ is a smoothed version of the spectrum $\vec{U}_{i}$.

In addition to this, the user can choose to use a further refinement to the unfolding method known as Compton subtraction [26]. It is used to further control the fluctuations in the unfolded spectrum. The basic concept behind it is to use the previously unfolded spectrum to decompose $\vec{R}$ into parts corresponding to the full-energy, single and double escape and annihilation peaks, and the ‘‘rest’’ which comes from Compton scattering and similar processes. Each of these parts, save for the full-energy peak, are then smoothed with the detector resolution before they are subtracted from $\vec{R}$. The resulting spectrum normalized to maintain the number of counts. The idea is that this should give an unfolded spectrum with the same statistical fluctuations as in the original spectrum $\vec{R}$.

Appendix B The first-generation method

In this appendix we describe the idea behind the first-generation method of Ref. [27] and its implementation in OMpy. Let $FG(E_{\gamma})_{E_{x}}$ denote the first-generation $\gamma$-ray spectrum, i.e., the intensity distribution of $\gamma$-ray decay from a given excitation energy $E_{x}$, as function of $\gamma$-ray energy $E_{\gamma}$. Generally, the nucleus will decay from $E_{x}$ down to the ground state by emitting a cascade of $\gamma$ rays, which forms the total $\gamma$-ray spectrum. The total, or all-generations $\gamma$-ray spectrum, denoted $AG(E_{\gamma})_{E_{x}}$, can be viewed as a superposition of the first-generation spectrum and a weighted sum of the all-generations spectra of excitation energies below,

[TABLE]

Here, $w(E_{x}^{\prime})_{E_{x}}$ is a weight factor that gives the decay probability from $E_{x}$ to $E_{x}^{\prime}$, and $n(E_{x}^{\prime})_{E_{x}}$ is a normalization factor which corrects for the varying cross section to populate the $E_{x}^{\prime}$ bins. Note that the normalization factor $n(E_{x}^{\prime})_{E_{x}}$ times the all-generations spectrum $AG(E_{\gamma})_{E_{x}^{\prime}}$ gives the spectra of $\gamma$ rays $\left\langle AG(E_{\gamma})_{E_{x}^{\prime}}\right\rangle$ emitted from the $E_{x}^{\prime}$ bin for each single population.

Normalization and multiplicity estimation

The normalization factor $n(E_{x}^{\prime})_{E_{x}}$ can be estimated from the singles spectra $S(E_{x})$, which contain the number of reactions populating the excitation energy bin $E_{x}$ and thus the number of $\gamma$-ray cascades out of this level,

[TABLE]

However, usually it is determined from the total $\gamma$-ray spectrum by the relation

[TABLE]

where $\langle M(E_{x})\rangle$ and $N(E_{x})$ denote the average $\gamma$-ray multiplicity and the total number of counts, respectively, at excitation energy $E_{x}$. This reformulation uses the fact that $S(E_{x})=N(E_{x})/M(E_{x})$. In OMpy there are two ways to determine the average multiplicity $\langle M(E_{x})\rangle$. The initial idea, the total multiplicity estimation is given in Ref. [9] as

[TABLE]

where $\langle E_{\gamma}\rangle$ is the weighted-average $\gamma$-ray energy at excitation energy $E_{x}$. Due to the detector threshold, we are usually not able to measure all $\gamma$ rays, and this will artificially increase $\langle E_{\gamma}\rangle$. To solve this problem, a statistical multiplicity estimation has been added to the oslo-method-software and is adapted in OMpy. The underlying idea is that in heavier nuclei, like in the rare-earth region, $\gamma$ rays from entry states at higher excitation energy will decay down to the yrast-line, where it enters at an energy denoted here as $E_{\mathrm{yrast}}$. From there, the $\gamma$ rays follow a non-statistical decay to the ground state. For heavy nuclei, there are many levels at low excitation energies, such that it is assumed here that the yrast transitions will proceed with many $\gamma$ rays of so low energy, that they are usually below the detector threshold. In that case, we can replace the excitation energy in Eq. (24) by the apparent excitation energy $\tilde{E}_{x}=E_{x}-E_{\mathrm{yrast}}$. A challenge in this method is to correctly estimate the entry energy $E_{\mathrm{yrast}}$. Thus we recommend to use the total multiplicity estimation whenever the experimental conditions allow for it.

Weight function and iteration

The weight function $w(E_{x}^{\prime})$ encodes the probability for the nucleus to decay from $E_{x}$ to $E_{x}^{\prime}$, and is in fact nothing but the normalized first-generation spectrum for $E_{x}$,

[TABLE]

By rewriting Eq. (B), we obtain

[TABLE]

This is a self-consistent set of equations for the $FG$ spectra, which we solve by an iterative procedure, starting with a set of trial functions $FG(E_{\gamma})_{E_{x}}$ and iterating until convergence is reached. In OMpy, the trial functions are chosen as constant functions, i.e.  with the same value for all $E_{\gamma}$. In the original implementation of the first-generation method, the trial functions are instead chosen based on a Fermi gas level density model [85, 61]. We have checked with OMpy that this gives identical results as with constant functions. The iterative procedure may sometimes lead to a negative number of counts in the resulting first generation spectra. As for the unfolding method, it is not clear whether they originate from the statistical nature of the data or from inaccuracies in the iterative procedure. The current default solution in OMpy is to remove negative bins from the obtained FG spectra. However, as described in Section 7, the user can choose to rather redistribute the negative counts to neighboring bins using the Matrix.fill_negative option. In cases where a few bins obtain large, negative counts, we advise the user to carefully analyze the results before processing them further.

Based on the development in oslo-method-software a small variation was added to OMpy to ensure better convergence for higher iteration numbers. Starting from iteration 5, each iteration $i$ of the $FG$ spectra is calculated at $FG_{i}=0.7{FG}^{\prime}_{i}+0.3FG_{i-1}$, where ${FG}^{\prime}_{i}$ is the spectrum from the analytical solution of (26). This is comparable to a fixed learning rate in machine learning algorithms.

Appendix C Derivation of the Oslo method equation

Here, we derive the relationship between the distribution of primary $\gamma$ rays, and the strength function and level density. The derivation is based on Weisskopf [86, p. 214--217], and Blatt and Weisskopf [87, p. 342, 649], but generalizes the equations to take into account angular momentum conservation. This corresponds to the Hauser-Feshbach theory of statistical reactions [88], but we apply several simplifications, like the Brink-Axel hypothesis, such that we can arrive at the Oslo method equation, Eq. (1).

We start with Bohr’s independence hypothesis [29] for the cross section $\sigma$ in a nuclear reaction $a+A\to\mathrm{CN}^{*}\to c+C$.555More precisely, we assume that our reaction, e.g. a $(d,p)$ reaction, leads to a compound nucleus $\mathrm{CN}^{*}$, which subsequently decays by $\gamma$-ray emission. Let us denote all quantum numbers of the entrance channel, the compound nucleus and the exit channel by $\alpha$, $i$ and $f$, respectively. The cross section is then given by [87, p. 342], [29]

[TABLE]

where the compound nucleus formation cross-section $\sigma_{\mathrm{CN}}$ depends only on the entrance channel $\alpha$. Following Bohr, the decay probability of the compound nucleus $G_{\mathrm{CN}}$ depends only on the branching ratios of the compound nucleus level $i$ to a specific channel $f$, but not on the entrance channel. In the Oslo method, we select only excitation energies $E_{x}$ below the particle separation threshold, so the compound nucleus can decay by $\gamma$-rays only. The decay probability $G_{\mathrm{CN}}$ of the state $i$ is then simply the $\gamma$-ray branching ratio to a specific final state $f$ with energy $E_{f}$ and the spin-parity $J^{\pi}_{f}$,

[TABLE]

where $\Gamma_{\gamma}(i\to f)$ is a partial, and $\Gamma_{\gamma}=\sum_{f}\Gamma_{\gamma}(i\to f)$ is the total radiative width for the level $i$.

The spectrum of the decay radiation $n(i,f)\mathrm{d}E_{f}$ from the compound nucleus level $i$ is then given by the summation (or integration) of $\sigma(\alpha,f)$ over an interval $dE_{f}$ (up to a constant due to the flux of $a$ and density of $A$ which would cancel out later):

[TABLE]

where $\rho_{\mathrm{avail}}$ is density of accessible final levels at $E_{f}$, $C_{\alpha\to i}$ is a constant that depends only on the entrance channel (which determines the compound nucleus state $i$), and we have replaced the average partial radiative width $\langle\Gamma_{\gamma}(i\to f)\rangle$ by $\gamma$-ray strength-function666To keep standard notation we will denote both the $\gamma$-ray strength function and the final level by $f$. It should be clear from the context what we refer to. $f_{XL}$.

The $\gamma$-ray strength function $f_{XL}$ for a given multipolarity $XL$ and for decays from an initial level $i$ to final level $f$ is defined as [89]

[TABLE]

where $\langle\cdots\rangle$ denotes an average over individual transitions in the vicinity of $E_{x}$ (in practice defined by the energy binning resolution). This can be simplified using the dominance of dipole radiation ($L=1$) and a common generalization of the Brink-Axel hypothesis [49, 50], where the strength function is assumed to be approximately independent of $E_{x}$, $J$ and $\pi$,

[TABLE]

where we define the total dipole strength function $f_{1}$ as the sum of the electric and magnetic dipole strength, $f_{E1}$ and $f_{M1}$, respectively, $f_{1}=f_{E1}+f_{M1}$.

If we again use the dominance of dipole radiation in Eq. (C) (which leads to $\Gamma_{\gamma}\approx\Gamma_{\gamma}^{(L=1)}$), and assume parity equilibration of the level density, i.e.  $\rho(E_{x},J,+)\approx\rho(E_{x},J,-)$ we can write

[TABLE]

where $\rho(E_{x},J_{f},\mathrm{eq})$ denotes the level density of one parity, the notation emphasizing the assumption of parity equilibration777The selection rules dictate that dipole radiation changes the angular momentum $J$ by at most one unit. For $M1$, the parity is unchanged, while for $E1$ it flips. This determines the density of available final levels for the decay. In the case of $J_{i}=1/2$ the sum runs over $J_{f}=\{1/2,3/2\}$, and in the case of $J_{i}=0$, the sum only runs over $J_{f}=1$, since $J=0\to J=0$ transitions are forbidden. and $C^{\prime}_{\alpha,i}=C_{\alpha,i}/\Gamma_{\gamma}$. We may write the partial level density $\rho(E_{x},J,\mathrm{eq})$ as

[TABLE]

where $g$ denotes the intrinsic spin distribution of the nucleus and $\rho(E_{x})=\sum_{J^{\pi}}\rho(E_{x},J,\pi)$ is the ‘‘summed’’ (or ‘‘total’’) nuclear level density. With this, we can further simplify the sum over the final levels in Eq. (C):

[TABLE]

which is a good approximation except for the case of $J_{i}=0$ (and $J_{i}=\nicefrac{{1}}{{2}}$), where the selection rules allow only transitions to $J_{f}=1$ ($J_{i}=\{\nicefrac{{1}}{{2}},\nicefrac{{3}}{{2}}\}$) states.

Next, we write $n(i,f)\mathrm{d}E_{f}$ more explicitly as $I(E_{i},J^{\pi}_{i},E_{\gamma})$ and exploit probability conservation,

[TABLE]

where $P_{J_{i},\pi_{i}}(E_{i},E_{\gamma})$ is the probability to decay from an initial excitation energy bin $E_{i}$ with a $\gamma$ ray of energy $E_{\gamma}$, the subscripts $J_{i}$ and $\pi_{i}$ limit the initial levels to of a given spin and parity, respectively, and $C_{E_{i},J^{\pi}_{i}}$ is a normalization constant. Note that the normalization constant $C^{\prime}_{\alpha,i}$ cancels out (which includes compound nucleus cross-section $\sigma_{\mathrm{CN}}$, the total radiative width $\Gamma_{\gamma}$, and the density of intrinsic levels $\rho(E_{i},J^{\pi}_{i})$).

The final step is to generalize this equation for the case where levels of different spins and parities are populated. Naively, we may just sum over the decays $P_{J^{\pi}_{i}}$ from all initial levels $J^{\pi}_{i}$

[TABLE]

from which we would already recover the standard Oslo method equation for a suitable normalization constant $C_{E_{i}}$. However, we have to note that the normalization constants $C_{E_{i},J^{\pi}_{i}}$ in Eq. (35) depend on the spin and parity, and cannot be factored out. As we will see, this can be solved under the assumption that the cross-section $\sigma_{\mathrm{CN}}$ to create the compound nucleus at different spin-parities $J^{\pi}_{i}$ is proportional to number of levels in the nucleus (i.e. it is not spin-selective, but proportional to intrinsic spin distribution),

[TABLE]

Using Eq. (C) to (34), the (cross-section weighted) sum over the decay spectra from all populated levels is

[TABLE]

In principle, we also have to average over the excitation energy bin $E_{x}$. However, as the level density $\rho$ and the total average radiative width $\Gamma_{\gamma}$ are assumed to vary only slowly with energy, this will not lead to any changes in the equation above. The normalized spectrum is given by

[TABLE]

for a normalization constant $C_{E_{x}}$ that depends only on the excitation energy.

In the case of a more spin-selective population of the compound nucleus, like in $\beta$-decay one receives a weighted sum of the level densities. If we denote the normalized population of the levels per for each excitation energy by $g_{\mathrm{pop}}(E_{x},J,\pi)=\sigma_{\mathrm{CN}}(E_{x},J,\pi)/\sum_{J^{\pi}}\sigma_{\mathrm{CN}}(E_{x},J,\pi)$, Eq. (C) can be generalized as

[TABLE]

In summary, we have shown that the Oslo method equation arises naturally from the Bohr’s independence hypothesis for the compound nucleus under the assumption of i) the dominance of dipole radiation, ii) parity equilibration of the level density and iii) a compound nucleus cross-section $\sigma_{\mathrm{CN}}$ that is proportional to the spin distribution of intrinsic levels.

Appendix D Calculation of $\langle\Gamma_{\gamma}\rangle$

In this appendix we derive Eq. (6.1) used to calculate the average total radiative width $\langle\Gamma_{\gamma}\rangle$ from the level density $\rho$ and $\gamma$-ray strength function $f$. The derivation is based on Eqs. (7.19-7.23) in Ref. [87] and p. 106 in Ref. [50]. It bases on the same arguments as C, but for ease of readability, we reiterate the main points.

The radiative width $\Gamma_{\gamma}$ denotes the probability for a state to decay by $\gamma$-ray emission. In the multipolar expansion it is written as

[TABLE]

where $EL$ and $ML$ denote the electric and magnetic components of the radiation of order $L$, respectively. To simplify the derivation, we will now assume the dominance of dipole radiation, such that

[TABLE]

We will continue with the derivations for the $E1$ radiation, but similar equations hold for $M1$. Analogously to Eq. (C), we can express $\Gamma_{E1}$ as the sum of the partial decay widths $\Gamma_{E1,i\rightarrow f}$ from the initial level $i$ to allowed final levels $f$,

[TABLE]

where the initial level is at the excitation energy $E_{x}$ and has the spin $J_{i}$ and parity $\pi_{i}$. The angular momentum and parity of the final states are given by the usual selection rules.

We now rewrite the partial widths $\Gamma_{E1,i\to f}$ in terms of the strength function $f$, using the definition of $f$, Eq. (30), and the generalized Brink-Axel hypothesis, Eq. (C),

[TABLE]

where $\rho(E_{x},J_{i},\pi_{i})$ is the density of the initial levels, and the energy difference to the final state(s) is given by $E_{\gamma}=E_{x}-E_{f}$.

Next, we replace the the sum in Eq. (43) by an integral, where we note that the number of partial widths $\langle\Gamma_{E1}(E_{\gamma},E_{x})\rangle$ is proportional to the level density at the final states

[TABLE]

such that the average total radiative width $\langle\Gamma_{\gamma,E1}(E_{x},J_{i},\pi_{i})\rangle$ is given by

[TABLE]

where the $E1$ operator requires a change of parity for the final state, here denoted as $\bar{\pi}_{i}$. The same selection rules as given in footnote 7 have been applied, but we do not use the approximation of Eq. (34), as many nuclei have a initial spin $J_{i}$ of 0 or $\nicefrac{{1}}{{2}}$. Further, we assume parity equilibration of the level density, i.e.  $\rho(E_{x},J,+)\approx\rho(E_{x},J,-)$, and express the level density $\rho(E_{x},J,\pi)$ through the spin-distribution $g(E_{x},J)$, Eq. (33),

[TABLE]

At this point, to obtain the expression for $M1$ radiation, one only needs to exchange the $E1$ strength function $f_{E1}$ by the $M1$ strength function $f_{M1}$. Using Eq. (42), we find

[TABLE]

This is equivalent to Eq. (2.11) in Ref. [89] under the assumptions listed above. One can determine $\langle\Gamma_{\gamma}(E_{x},J_{i},\pi_{i})\rangle$ from neutron capture experiments, where usually slow neutrons are used, thus $E_{x}\approx S_{n}$. This is shown e.g. in Ref. [90], but we will repeat the derivation to get a comprehensive picture. The intrinsic spin of the neutron is $\nicefrac{{1}}{{2}}$, so with capture of order $\ell$ on a target (denoted with the subscript $t$) the entry states in the residual nucleus have the possible spins $J_{i}=J_{t}\pm\nicefrac{{1}}{{2}}\pm\ell$ and parity $\pi_{i}=\pi_{t}(-1)^{\ell}$. For s-wave capture on a target nucleus with $J_{t}=0$, there is only one possible $J_{i}$, and we can directly compare the experimental measurements to the calculations using (48). For $J_{t}>0$, often only the average over all resonance of different $J_{i}$ is reported. Using the level density $\rho(S_{n},J_{i})$ of the accessible $J_{i}$’s, we find

[TABLE]

where we substituted $\sum_{J_{i}}\rho(S_{n},J_{i})$ by the average neutron resonance spacing $D_{l}$ and used the assumption of the dominance of the dipole decay to rewrite the sum of the average $M1$ and $E1$ strengths as $f(E_{\gamma}$). Note that transitions for the highest $\gamma$-ray energies go to discrete states, not a quasi(-continuum). Thus, this integral is an approximation, and more precise calculations could distinguish between a sum over transitions to discrete states and the integral for the (quasi-)continuum region.