Study of a data analysis method for the angle resolving silicon telescope

A new data analysis method is developed for the angle resolving silicon telescope introduced at the neutron time of flight facility n_TOF at CERN. The telescope has already been used in measurements of several neutron induced reactions with charged particles in the exit channel. The development of a highly detailed method is necessitated by the latest joint measurement of the $^{12}$C($n,p$) and $^{12}$C($n,d$) reactions from n_TOF. The reliable analysis of these data must account for the challenging nature of the involved reactions, as they are affected by the multiple excited states in the daughter nuclei and characterized by the anisotropic angular distributions of the reaction products. The unabridged analysis procedure aims at the separate reconstruction of all relevant reaction parameters - the absolute cross section, the branching ratios and the angular distributions - from the integral number of the coincidental counts detected by the separate pairs of silicon strips. This procedure is tested under the specific conditions relevant for the $^{12}$C($n,p$) and $^{12}$C($n,d$) measurements from n_TOF, in order to assess its direct applicability to these experimental data. Based on the reached conclusions, the original method is adapted to a particular level of uncertainties in the input data.


Introduction
The neutron time of flight facility n_TOF at CERN is a neutron production facility aiming at measuring the neutron induced reactions. A massive lead spallation target irradiated by the 20 GeV proton beam from the CERN Proton Synchrotron serves as the primary source of neutrons, delivering an extremely luminous white neutron beam spanning 12 orders of magnitude in energy -from 10 meV to 10 GeV. The n_TOF facility features two experimental areas: Experimental Area 1 (EAR1), horizontally placed at 185 m from the spallation target, and the Experimental Area 2 (EAR2) vertically placed at 20 m above the target. While EAR1 is best adjusted to the high neutron energy and the high resolution measurements, EAR2 excels at the measurements with small, highly radioactive samples characterized by low cross sections for the investigated reactions. More details on the general features of the n_TOF facility and EAR1 itself may be found in ref. [1], while the specifics on EAR2 are addressed in Refs. [2][3][4]. An overview of the experimental program at n_TOF may be found in ref. [5]. A general overview of many different types of detectors used at n_TOF for the measurements of various types of the neutron induced reactions, together with the detailed description of the procedures for the analysis of electronic signals from these detectors, can be found in ref. [6].
A new, highly sophisticated silicon telescope (SITE) has recently been introduced at n_TOF for measurements of the neutron induced reactions with charged particles in the exit channel [7]. It consists of two separate and segmented layers of 16 silicon strips, 5 cm × 3 mm each, placed in these two telescopes as front and rear, the front one being parallel to the sample and covering the forward angles, with the rear one being parallel to the neutron beam and mostly covering the backward angles (see figure 1b). The analysis of the experimental data on the 12 C(n, p) and the 12 C(n, d) reaction is under way, pending the development of a new analysis method for extracting the relevant reaction parameters. Most important among these is the absolute cross section. The angle-differential cross sections for the reaction flow via the separate excited states of a daughter nucleus ( 11 B [14] or 12 B [15]) would also be highly desirable. However, the reliable decoupling of these states might not be possible at the level of statistics expected from the latest measurement.
The purpose of this paper is twofold. The first is to develop the necessary formalism for the analysis of the data obtained with the multi-channel telescope (section 2). The second is to investigate its applicability on the artificially generated dataset resembling the first experimental dataset from n_TOF to which the the method is to be applied at a later date: that of the 12 C(n, p) reaction (section 3). In doing this we aim (1) to provide the future users of the method with all the necessary steps and considerations to be taken into account in extracting the optimal set of physical parameters from a given measurement; (2) to provide an honest assessment of the direct applicability of the method to a dataset of a given level of uncertainties, in particular the one expected from the 12 C(n, p) measurement, and (3) to provide alternative solutions in case the direct application proves to be unreliable due to the level of uncertainties in the extracted results (section 4). Since we develop the method having a specific 12 C(n, p) reaction in mind, it cannot be overemphasized that the procedure is aimed at and designed for the particular detector setup, based on the ∆E-E telescoping principle, rather than for the particular reaction of even the type of reaction. Therefore, at no point should the method be considered as limited to this specific (type of) reaction, nor should the conclusions regarding a particular 12 C(n, p) measurement from n_TOF be misinterpreted for some general limitation of the method itself.

Method derivation
Let θ be the angle of proton emission in the center of mass frame (of the incoming neutron and 12 C nucleus before the reaction, and of the outgoing proton and the 12 B nucleus after the reaction), relative to the direction of the neutron beam. We immediately introduce: as a relevant variable. For simplicity of terminology we still refer to χ as the angle of the proton emission. Let N i j be the total number of protons detected in coincidence by the (i, j)-pair of strips, with the first index i denoting some of the thin ∆E-strips and the a second index j denoting some of the thick E-strips from any telescope (front or rear). Let E be the energy of the incident neutron. The proton produced by the neutron of sufficiently high energy might be emitted leaving the 12 B nucleus in any of the energetically accessible states. Thus, the proton energy is clearly contingent on the daughter nucleus' excited states. Denoting these states by x (x = 0 being the ground state), we define the probability ε i j (x, E, χ) for the coincidental detection -by the (i, j)-pair of stripsof protons produced by the neutrons of energy E and emitted under an angle χ leaving the 12 B nucleus in a state x: with d 2 N i j (x, E, χ) as the number of the detected protons and d 2 N (x, E, χ) as the number of protons emitted under such conditions. These probabilities may easily be obtained from the dedicated simulations, described in appendix A. It should be noted that they reflect the properties of the experimental setup itself, and are independent of the angular distribution of the emitted protons.
Only for illustration purposes, figure 2 shows the coincidental detection probabilities ε i j (0, 20 MeV, χ) for protons produced by 20 MeV neutrons, leaving the 12 B nucleus in the ground state (x = 0). The probabilities are shown for some arbitrarily selected, i-th ∆E-strip in coincidence with the several closest E-strips. During the method implementation these curves, i.e. their smooth(ed) forms never have to be constructed, as their integrals can be calculated as the weighted sum of the simulated counts. The issue is further addressed in appendix A.
The number of protons emitted under the described specific conditions may be decomposed as: with φ(E) as a time-integrated energy dependent neutron flux irradiating the sample: φ(E) = dΦ(E)/dE, dΦ(E) being the total number of neutrons of energy E intercepted by the sample. The multiple scattering factor µ(E) describes an increase in the neutron flux at an energy E due to the energy loss of higher-energy neutrons by means of the multiple scattering inside the sample itself. With (x, E, χ) as the partial cross section for the 12 C(n, p) reaction, i.e. for a particular reaction of interest, Σ tot (E) is the total cross section for any neutron induced reaction in the carbon sample. Finally, η is the areal density of the sample, as encountered by the neutron beam, in the number of atoms per unit area. While the term 1 − e −ηΣ tot (E) gives a probability for any neutron reaction to occur (the exponential term itself being the transmission probability), the differential ratio (x, E, χ)/Σ tot (E) governs the probability of that reaction being the one of interest.
The differential cross section may now be decomposed as: with σ(E) as the total cross section for the 12 C(n, p) reaction, ρ(x, E) as the energy-dependent branching ratios for the reaction flow via the particular excited state of 12 B, and A(x, E, χ) as the angular distribution of protons specific to that state. From eq. (2.3) we now isolate all the terms that are independent of the detector setup, while being available from the experiment, simulation or any evaluation database: The neutron flux φ(E) at EAR1 (as well as the flux at EAR2 [16]) is available from the dedicated measurements at n_TOF [17]. Even in a general case of a thick sample, the multiple scattering factor could be obtained from the dedicated simulations if the total cross section Σ tot (E) and the elastic scattering cross section Σ el (E) for carbon were known with sufficient precision, which they are [12]. However, as the very thin carbon sample was used during the energy-differential measurement -0.25 mm [13], with the thickness of 0.35 mm, i.e. an areal density of η = 4 × 10 −3 atoms/barn being intercepted by the neutron beam due to the sample tilt of 45 • (figure 1b) -a thin sample approximation becomes highly appropriate. In this approximation the deviation of the multiple scattering factor from unity is completely negligible: µ(E) ≈ 1, while the full fractional term from eq. (2.5) approximates to η due to ηΣ tot (E) 1 within the entire neutron energy range of interest. Hence: w(E) ≈ ηφ(E). 3) may now be rewritten as: We now take into account that due to the energy spread of the neutron beam the experimental data must be analyzed within the energy intervals of finite width. We use the following notation for one such interval: meaning that all the later quantities denoted by E are either integrals or averages over E, or that they may be separately and independently selected for each such energy interval. Since any particular method implementation requires a weighted averaging over w(E), we immediately introduce the following norm: Returning to the differential number of protons d 2 N i j (x, E, χ) detected by a particular pair of strips and recalling that there may be multiple excited states of the daughter nucleus contributing to the reaction, we may write the expression for the total number of protons detected by the (i, j)-pair of strips: with X E as the highest excited state affecting the data from the energy interval E. It should be noted that the total detected counts N (E) i j taken for analysis will also be dependent on the energy deposition cuts imposed on the experimental data. We will consider this dependence implicitly absorbed within the terms N (E) i j and, consequently, the corresponding detection probabilities ε i j (x, E, χ). We now define an arbitrary bijective mapping: allowing us to write eq. (2.10) in a single index, which will soon become useful in bringing the system to an appropriate matrix form. This bijection never needs to be explicitly constructed. Using this formal manipulation in conjunction with applying Eqs. (2.2) and (2.7) to eq. (2.10), we arrive at the master equation: Our goal is now to bring this equation into the matrix form: by constructing the vector ì N (E) of total detected counts N (E) α , by designing an appropriate matrix E E and by isolating the sought parameters of the partial cross section within the vector ì P (E) . We will obtain this matrix form by decomposing the angular distributions into partial waves (Legendre polynomials).
Before proceeding further let us put forth the tools and considerations common to any particular implementation of the method. Let us denote by R E the number of relevant pairs of strips composing the experimental dataset from ì N (E) and by P E the number of the partial cross section parameters from ì P (E) . Then we can select at most R E parameters to reconstruct: (2.14) When P E < R E , the best solution to this system may be found by means of the weighted least squares method [18]: with V E are the covariance matrix of the input data, allowing for the propagation of experimental uncertainties in order to obtain the covariance matrix V E of the reconstructed parameters and their respective uncertainties δP (E) β as: From the raw results obtained from eq. (2.15) we will have to calculate various consequent quantities -the total reaction cross section, branching ratios and the angular distribution parameters -while dealing with the high uncertainties and possible correlations in the reconstructed ì P (E) . Therefore, we are well advised to take into account the effects of the full covariance matrix upon the propagation of uncertainties. Let any of these quantities be the scalar function of ì P (E) that we generally denote as: F E ≡ F ì P (E) . Then the respective uncertainty δF E may be expressed as: with J F indicating the conventionally defined Jacobian matrix of the function F.

Partial waves decomposition
We start by decomposing the angular distributions into the selected number of partial waves, i.e. Legendre polynomials P l ( χ): with the maximum wave number L E (x) freely adjustable to any given excited state, within the constraints of the total number R E of the available data points. Eq. (2.12) may now be rewritten as: where we recognize the appearance of the weighted averages · E and · l , with w(E) and P l ( χ) as the respective weighting functions. We now approximate the average product by the product of averages: with the single and double averages appearing as: In analogy to eq. (2.11) we introduce another arbitrary bijective mapping: never having to be explicitly constructed, but allowing for a unique-index labeling. In that, β spans the range of all free parameters, i.e. the total number of coefficientsā (E) xl : β = 1, ..., P E . As it holds: , from eq. (2.14) we have the following constraint upon the distribution of Legendre coefficients among the relevant exited states: ) is now recast as: having thus been brought into the matrix form from eq. (2.13), with the appropriate definitions: The entire solution ì P (E) and the corresponding uncertainties are now easily found from Eqs. (2.15) and (2.16).
Applying the normalization condition i.e. all the 0 th terms are fixed and carry the entire angular distribution norm. Plugging this result into eq. (2.29): x /2 and combining it with the normalization condition X E x=0ρ (E) The next step consists of identifying the branching ratios as: culminating in the separation of the angular coefficients: The uncertainties δσ (E) , δρ (E) x and δā (E) xl follow directly from eq. (2.17), according to the full covariance matrix V E from eq. (2.16) .
When the total number of the relevant excited states becomes so large that the total number P E of required parameters becomes comparable to the total number R E of available data points, and/or when these points are affected by large uncertainties, the coefficientsā (E) xl exhibit substantial uncertainties themselves and the contributions from the particular excited states can not be reliably separated. In this case one may attempt to reconstruct the "global" partial wave coefficients averaged over the excited states:ā hoping for a manageable uncertainty in the total contribution to a given partial wave. The factors f ,defined as: simply take into account whether a given partial wave was adopted for a given excited state. As all the 0 th terms are fixed by eq. (2.30), we immediately haveā (E) 0 = 1/2 and δā (E) 0 = 0.

Method implementation
We illustrate the implementation of the method by applying it to the 12 C(n, p) data artificially generated by the Geant4 simulations. We consider here the data from 1 MeV wide energy range , approximately where this reaction's cross section is expected to reach its maximum. The branching ratios and the angular distributions for each relevant excited state were arbitrarily constructed.

Selecting the excited states
The excited states contributing to the reaction within the given neutron energy range E need to be clearly identified, as the method requires them to be known in advance. For the 12 C(n, p) reaction, there are total of 15 states in the 12 B daughter nucleus with the energy threshold E thr below the upper limit of the considered neutron energy range (E thr < 20.5 MeV) [15]. Their excitation energies together with the corresponding Q-values and the energy thresholds in the laboratory frame are listed in table 1. While all these states contribute to the reaction, not all of them necessarily contribute to the totality of the detected counts, especially those very close to the reaction threshold. The reason is threefold: (1) the very low reaction cross section close to the threshold; (2) the pronounced forward boost of the produced protons in the laboratory frame, making them miss most of the detection setup; (3) the low proton production energy causing them to be stopped by the sample itself, never reaching the detectors at all. Therefore, it needs to be estimated in advance which states may be excluded from the analysis of the experimental data. As the exact evaluation of the expected amount of the detected counts from each state is, of course, impossible without the prior knowledge of the partial cross sections for each separate state (their branching ratios and angular distributions), one needs to rely on some reasonable estimate. One such useful figure of merit is the approximate probability estimatorε (E) x for the coincidental detection by any pair of the ∆E-E strips:ε constructed by assuming -in the absence of any prior information -the isotropic angular distribution of protons in the center of mass frame: A(x, E, χ) ≈ 1/2, and applying the same energy  [15], the corresponding Q-values and the energy thresholds E thr for the 12 C(n, p) reaction in the laboratory frame.
x deposition cuts as are to be applied to the experimental data. Figure 3 shows thus obtained probability estimates for the relevant 12 B states. Although the portion N (E) x of the produced protons still remains entirely unknown, the observed decrease inε (E) x together with the expected decrease in N (E) x for the higher states allows one to make informed estimates about the relevance of the expected partial contributions N (E) x to the detected counts: x . From these considerations applied to figure 3 we elect to include only the first 11 states (up to the 10 th excited one, i.e. X E = 10) for further analysis. The artificial data to be analyzed were, of course, simulated by including all 15 states with the energy thresholds below the upper limit of the considered neutron energy range.
It must be pointed out that this exclusion of higher states from the analysis may, in principle, affect the cross section normalization, as the branching ratios of the excluded states become unobtainable. However, as already discussed, the cross sections around the energy thresholds for these states are expected to be negligible and so is their impact upon the total reaction cross section. Still, if there were reasonable indications to the contrary, one should be aware that the reconstructed cross sectionσ (E) is only partially contributed by those states that were kept for the analysis.

Varying the amount of partial waves
The highest wave numbers L E (x) for each excited state are evidently the method's adjustable parameters. For the total of R E relevant pairs of strips from eq. (2.14), there is a total of R E X E +1 selections of L E (x) satisfying the constraint from eq. (2.26), with · · denoting a binomial factor. For R E = 60, as used later, and X E = 10 this amounts to approximately 3.4 × 10 11 combinations. If we were to impose some maximum admissible wave number L E that may be assigned to any particular state -implying, for purpose of these simple estimates, that the selection of L E itself must be such that (L E + 1)(X E + 1) ≤ R E , in order for each of X E + 1 states to be allowed L E + 1 waves -then the number of possible selections for L E (x) reduces to (L E + 1) X E +1 . For example, the maximum value L E = 4 compatible with R E = 60 and X E = 10 yields approximately 4.9 × 10 7 combinations. However, the following physical argument helps us in reducing the number of possible combinations even further, by keeping only the physically sensible selections of L E (x). We consider that close to the reaction threshold the nuclear reactions are expected to be isotropic (in the center-of-mass frame), while the anisotropy is expected to appear (and possibly intensify) with increasing incident particle energy. This suggests that the higher excited states -characterized by a higher threshold -should not be assigned more partial waves than the lower states, i.e.: For the maximum admissible wave number L E , the number of combinations consistent with this constraint is now reduced to L E +X E +1 X E +1 . For example, the maximal value L E = 4 compatible with R E = 60 and X E = 10 leaves the total of 1365 combinations. All we need now is the algorithm for constructing such combinations. For the maximum wave number L E to be assigned to any state, the particular combination of nonincreasing L E (x) values may be uniquely defined by the set of L E states Λ ( = 1, . . . , L E ) at which the maximum wave number L E (x) increases by 1. In other words, Λ form a set of states such that L E (x) = ends at x = Λ , i.e.: with additional fixed boundaries Λ 0 = X E and Λ L E +1 = −1, defined for the convenience of the implementation. The algorithm now reduces to generating all combinations (L E -tuples) of Λ such that: It is easy to confirm that if Λ = −1 for all , then L E (x) = 0 for all x, i.e. all the states are assigned only the isotropic component. The other extreme is Λ = X E for all , meaning that L E (x) = L E for all x, i.e. all the states are assigned the maximum allowed number of partial waves.

Optimizing the model parameters
The obvious question now is how to select an optimal combination of the wave numbers L E (x).
We propose here a simple -and by no means unique -selection principle. As the variations in L E (x) directly affect the number of the model parameters: , the reduced chi-squared estimator X 2 lends itself easily to a quick and efficient evaluation of the goodness of the fit: The rightmost expression holds when the covariance matrix V E of the input data is diagonal, i.e. when the correlations between the components of ì N (E) are negligible. As opposed to the goodness of fit -which will for large R E systematically improve by increasing the number of partial waves, as long as P E does not closely approach R E -the reliability of the fit, reflected through the uncertainties in the reconstructed ì P (E) , rapidly degrades with increasing number of parameters. For estimating this reliability we propose a simple calculation of the uncertainty δX 2 in the chi-squared value from eq. (3.5) by means of eq. (2.17), since X 2 is sensitive to all the fitted parametersunlike, for example, the reconstructed cross sectionσ (E) from eq. (2.31). In the context of our problem the minimization of X 2 and its uncertainty δX 2 seem to be opposing objectives. Therefore, we propose to minimize their product X 2 δX 2 as the simplest estimator that should at its minimum provide the optimal tradeof between the goodness and the reliability of the fit.
There are additional issues to consider. For the number of partial waves too inadequate for a given set of the experimental data, some of the branching ratiosρ (E) x from eq. (2.32) may turn out to be negative or greater than unity -a clear signature of the badness of the fit, going beyond the particular values of X 2 . These fits should be immediately rejected as physically unsound, i.e. disqualified from any kind of optimization procedure, be it the minimization of X 2 δX 2 or some alternate technique.
Yet another quality control mechanism consists of checking if the reconstructed angular distributions for each angular state: become negative at any point. If so, such fits may also be immediately rejected, regardless of their goodness. One should be wary, however, in making such rejections when there are prior indications that some states may indeed feature the very low branching ratios or highly anisotropic angular distributions that locally come close to zero. In this case any statistical fluctuation in the input data may easily drive the reconstructed results toward the negative values, while the results do remain reasonably reliable representations of the true reaction parameters. It should be noted that the reconstructed branching ratios may discard all fits as unphysical, if every combination of wave numbers L E (x) produces at least one negativeρ (E) x . On the other hand, the isotropic angular distributions will always pass the negativity test, so that the fully isotropic fit (L E (x) = 0 for all x) will always be accepted, based on the positivity of the angular distributions themselves.
Prior to calculating X 2 , δX 2 or any consequent quantity to be used in judging the suitability of the fitted result, one may also consider manually eliminating from the set of fitted parameters those P (E) xl that, according to eq. (2.33), yield the angular coefficients too small (|ā (E) xl | 1) or unreasonably large (|ā (E) xl | 1) in magnitude. For the associated β, this is most easily done by setting P (E) β = 0 and [V E ] ββ = [V E ] β β = 0 for all β within the covariance matrix from eq. (2.16). This procedure helps in regularizing the fit, as the exceedingly small |ā (E) xl | are commonly the sporadic results caused by the finite precision data, while the distinctly large |ā (E) xl | are expected to appear as the consequence of overfitting the statistical fluctuations in the input data. One should, of course, be prepared for the closer inspection and the critical evaluation of the results if the optimal set of parameters happens to be precisely thus manipulated set. However, what is expected from this procedure is the artificially induced increase in the fit suitability estimator X 2 δX 2 , such that some alternative set of parameters takes precedence as the optimal one.
In summary, we propose to identify the optimal combination of the maximum wave numbers L E (x) by minimizing the product X 2 δX 2 -or any such estimator balancing between the goodness and the reliability of the fit -while taking into account the physical soundness of the results, whether by immediately rejecting those physically inadmissible or by appropriately penalizing them during the optimization procedure.

Method investigation: 12 C(n,p) data
We now test the method on a particularly challenging example of the artificially generated 12 C(n, p) data, as means of appraising its applicability to the experimental data from n_TOF. The simulated dataset -the set of counts N (E) α detected by a particular pair of strips -was obtained from the same Geant4 simulations as used for obtaining the coincidental detection probabilities, i.e. the central design matrix E E . The neutron energies were sampled within the 1 MeV wide interval E = [19.5 MeV, 20.5 MeV], all 15 states from table 1 were used in constructing the dataset, while only the first 11 states from figure 3 were considered for the reconstruction. For the buildup of the test counts an arbitrarily constructed branching ratios ρ(x, E) for each of the 15 states were used (represented by later figure 5), together with the angular distributions A(x, E, χ) arbitrarily constructed for each state, which were all designed from the three lowest Legendre polynomials (P 0 , P 1 , P 2 ). Figure 4 shows the relevant set of coincidental counts recorded by different ∆E-E pairs of strips, ordered by magnitude. While there are (16 E-strips)×(16 ∆E-strips )×(2 telescopes) = 512 possible pairs of strips in the used SITE configuration from figure 1b, one can see from figure 4 that only the tenth of those are characterized by a sufficient coincidental detection probability to be considered for analysis. It should be noted that the counts from figure 4 were constructed from an exceedingly large dataset, featuring the negligible statistical fluctuations. In order to easily generate the statistical variations in the dataset to be taken for analysis, we first scale these counts to a desired level of statistics (thus constructing their statistically expected values) and then generate the appropriate Poissonian fluctuations. For purposes of this demonstration we keep only those coincidental counts N (E) α that are higher than 5% of the maximum value (the dashed threshold from figure 4). Depending on a particular realization of the Poissonian fluctuations, around R E = 60 relevant pairs of strips meet this condition. The statistical uncertainties are then assigned to these counts by setting the diagonal elements of the input covariance matrix V E from Eqs. (2.15) and (2.16) In order to vary the maximum wave numbers L E (x) for each excited state we follow the procedure from section 3.2, adopting the maximum supported value L E = 4. We choose the number of counts from the most efficient pair of strips to be: max N (E) α = 10 6 , making the total number of counts detected across all kept pairs: R E α=1 N (E) α = 2 × 10 7 . The reason behind this selection is rather simple and carries the critical repercussions for the analysis of the experimental data from n_TOF: at lower statistics basically all the fits are discarded due to the appearance of the negative branching ratios. In other words, for almost all generated instances of Poissonian fluctuations all the fits (for any combination of state boundaries Λ ) produce at least one negativeρ (E) x . One must be careful at this point not to confuse max N (E) α = 10 6 with some minimum intrinsic level of reliable statistics. Instead, it reflects an amount of excited states at play: a high number of states naturally requires high statistics if they were to be reliably disentangled one from the other.
We now appraise the method based on the accuracy and uncertainty of the reconstructed parameters. For a condensed demonstration of the results on the reconstructed angular distributions,  Figure 5: Typical example of the reconstructed set of branching ratios, recovered by an optimal set of wave numbers assigned to each excited state. Only the first 11 states were considered for the reconstruction, as the rest of them hardly contribute to the detected counts or not at all.
we we use the overall distribution A E ( χ), averaged over all excited states: The reference distribution stems from the arbitrarily constructed distributions A(x, E, χ) for each of the 15 states contributing to the reaction (X + E = 14; see table 1). The reconstructed distribution, as denoted by , is contributed by the reduced number of states taken for the analysis (X E = 10). After applying the method to the different realizations of Poissonian fluctuations, our conclusions are rather simple and straightforward. The fits yielding an unphysical set of branching ratios also grossly misidentify the overall angular distribution A E ( χ) and, in general, the reconstructed cross sectionσ (E) , reflecting the absolute normalization of the data. As such, they should indeed be immediately rejected. Among the physically admissible fits (if there are any at all) the ones identified as optimal do seem to reasonably reconstruct both the overall angular distribution and the cross section, at least under the level of statistics adopted here out of necessity. However, the set of reconstructed branching ratios themselves most often seems to be unrepresentative of the true results, as illustrated by a typical example from figure 5. The example from figure 5 also shows that their uncertainties may also be grossly underestimated and unrepresentative of their error. Therefore, the reconstructed branching ratios should be taken with maximum caution.
At the adopted level of statistics most often there seems to be little difference in the results obtained by minimizing X 2 or the proposed product X 2 δX 2 , as the physicality of the branching ratios serves as the main discriminator of the unreliable fits. Figure 6 shows an example when the difference in the reconstructed overall angular distributions obtained by minimizing X 2 or X 2 δX 2  Figure 6: Overall angular distribution recovered from an optimal set of wave numbers for each excited state, obtained by minimizing either X 2 (the goodness of the fit) or X 2 δX 2 (the tradeoff between the goodness and the reliability). turns out to be appreciable. This example clearly illustrates the superiority of optimizing the tradeoff between he goodness and the reliability of the fit. The power of this procedure lies not in reducing the uncertainties per se, but in penalizing the overfitting, i.e. in rejecting the sporadic parameters that unnecessarily and disproportionately increase the uncertainties in all other parameters, besides introducing their own excessive ones. Indeed, while the reference angular distribution from figure 6 was constructed as a linear combination of the 3 lowest Legendre polynomials, the one identified by minimizing X 2 allows for 5 of them (the maximum amount supported by L E = 4; a clear symptom of overfitting), while the minimization of X 2 δX 2 finds the combination of 4 partial waves as the optimal one.
Let us recall that with so many exited states at play, the physicality of the branching ratios serves as the primary discriminator of unreliable fits. For a significantly reduced number of states, this method of assessment becomes much more insensitive or even entirely unavailable in case of a single relevant, ground state. In that case the quality tradeoff between the goodness and the reliability of the fit remains the crucial, if not the only available method for identifying the optimal set of the fit parameters. Finally, at the adopted level of statistics the relative uncertainty in the reconstructed cross section σ (E) appears to be around 10%. As the statistically expected uncertainty scales as N (E) tot −1/2 with the total number N (E) tot of the detected counts, one can easily estimate the expected level of uncertainty at any level of statistics, provided that the available data produce any acceptable fit in the first place. Considering that the experimental n_TOF data on the 12 C(n, p) reaction are expected to provide 4 to 5 orders of magnitude less statistics than adopted for this demonstration [13], even if they could be fitted without all fits failing the physicality test, the uncertainty inσ (E) is thus expected to be at least an order of magnitude grater than reconstructed cross section itself! Hence, the direct application of the full reconstruction method presented up to this point is ill-adjusted to these experimental data, due to the particularly unfavorable combination of the available statistics and the amount of excited states at play. This outcome should not be confused with some intrinsic shortcoming of the method itself, as there is a limit to the quality of the results that could be extracted from the data of a finite statistical precision. Fortunately, this eventuality was foreseen in advance of the experiment and the experimental setup was specially optimized so as to minimize the systematic effects due to the alternative approach to the analysis of these data. This approach consists of utilizing the reduced variant of the method, by adopting a priori information on the branching ratios and the angular distributions from an outside source -such as the TALYS theoretical model [19], adjusted to the preexisting experimental data -and aiming solely at the reconstruction of the absolute cross sectionσ (E) . This reduced variant is addressed in the following section.

Method reduction
Even when the full unfolding procedure may not be meaningfully applied due to the uncertainties in the input data limiting the usefulness of the output results, the method formalism from section 2 still remains relevant, as it clearly establishes the connection between the measured observables and the underlying reaction parameters. Furthermore, the coincidental detection probability of the experimental setup must be characterized -most appropriately by means of the dedicated simulations described in appendix A -regardless of the particular approach to the data analysis. Starting from eq. (2.12), one may derive any particular variant of the unfolding procedure, be it the reduction of the one from section 2.1 or even some further extension, shortly addressed in appendix B. Motivated by the status of the experimental n_TOF data on the 12 C(n, p) reaction, we consider here the adoption of a priori information on the branching ratios and angular distributions, aiming solely at the reconstruction of the absolute cross section. Assuming that information to be available from an independent source, eq. (2.12) may be linearized as: with the vector ì (E) (as a matrix of a reduced dimensionality) standing in place of the design matrix E E from eq. (2.13) and the single unknownσ (E) replacing the previous set ì P (E) of underlying reaction parameters. While the definition ofσ (E) stays the same as in eq. (2.21), ì (E) is now defined by components as: where the branching ratios ρ(x, E) and angular distributions A(x, E, χ) are taken from an outside source of information. Applying Eqs. (2.15) and (2.16) -while taking the covariance matrix V E to be diagonal and composed of the uncertainties δN (E) α in the detected counts: -the final solution for the sought cross section may now be written in a rather simple closed form: with the associated uncertainty: It should be noted that this procedure still makes full use of all the experimentally available information from separate pairs of ∆E-E strips. This feature is in clear opposition with the more extreme reduction of the method, taking only the total number N (E) tot = R E α=1 N (E) α of coincidental counts detected across an entire detection setup, in conjunction with its total detection probability (E) α in order to obtain the absolute cross section overly simplistically as tot / (E) tot , thus defeating any benefit of having used a high-end telescope -in particular, its sophisticated dissociation into multiple strips.
Evidently, the main challenge with thus reduced method is the estimation of the systematic uncertainties brought on by the out-of-necessity adopted branching ratios and angular distributions. An indication of those uncertainties -and a conservative one, at that -may be obtained by adopting all the involved angular distributions as isotropic: A(x, E, χ) = 1/2, and recalculatinḡ σ (E) . The difference between the externally provided and all-isotropic distributions is to be taken as representing the extreme case of the possible disparity with the true angular distributions. Another possibility is taking among the externally provided distributions only the branching ratios or the angular distributions as given, and unfolding the data with the other type of distributions unconstrained. Comparing these alternative results forσ (E) allows for an informed estimate of the systematic uncertainties.

Conclusions
A new angle resolving stripped silicon telescope (SITE) has recently been introduced at the neutron time of flight facility n_TOF at CERN for the measurements of the neutron induced reactions with the charged particles in the exit channel. Its outstanding detection properties have already been demonstrated in the challenging measurement of the 7 Be(n, p) reaction, relevant for the famous Cosmological Lithium Problem. The joint energy-differential measurement of the 12 C(n, p) and 12 C(n, d) reactions has also been recently performed at n_TOF, using the upgraded and specially optimized detector configuration consisting of the two separate silicon telescopes. As the nature of these reactions poses significant challenges for the meaningful data analysis -being affected by the multiple excited states in the daughter nuclei and featuring the anisotropic angular distributions of the reaction products -we have established a clear and detailed formalism behind the measured observables: the total number of the coincidental counts detected by any combination of ∆E-E pairs of silicon strips. From this formal connection we have developed and tested the unfolding procedure for the reconstruction of the underlying reaction parameters, consisting of the absolute reaction cross section, the branching ratios and the angular distributions of the reaction products for each excited state in the daughter nucleus. We have also addressed the finer points of the method implementation, thus providing the consistent and reliable methodology for obtaining the optimal set of the output parameters. Though the method may, in principle, reconstruct all these quantities separately, its performance may be severely limited by the amount of parameters -determined by the number of excited states and the level of anisotropy -as well as the level of uncertainties in the input data. By testing the method on the artificially generated dataset resembling the n_TOF measurement the of the 12 C(n, p) reaction, we have found little hope that the full unfolding procedure could be meaningfully applied to these particular experimental data, precisely due to the highly unfavorable combination of the large number of the excited states and the reduced level of statistics expected from the experiment. This unfortunate outcome should not be misinterpreted for the inherent deficiency of the method itself, as at some point all considered reaction parameters must be fully taken into account if the experimental data are to be properly described and reliably analyzed. Precisely the clarity of the formalism behind the method allows for its many alternative variants to be developed. One of these, to be applied to the measured 12 C(n, p) and 12 C(n, d) data, is the reduced procedure relying on the independent source of information on the branching ratios and angular distributions, aiming at the reconstruction of the absolute cross section as the central reaction parameter of interest. It is worth noting that thus reduced method still takes advantage of the distribution of the detected counts across the separate ∆E-E pairs of strips, as opposed to considering only the total number of counts across all of them. Thus retained angular sensitivity opens the possibility for the estimation of the systematic effects due to the adopted outside information (branching ratios and/or angular distributions), allowing for the informed assessment of the systematic uncertainties in the final results.

A Detection probability simulations
We describe the detection probability simulations and the use of the simulated data in the construction of the design matrix E E from eq. (2.28). For specificity, we again keep the 12 C(n, p) reaction in mind. The reaction details -except for its basic kinematics -are assumed unknown. The reaction itself or its basic details may not even be (properly) implemented in the used simulation package. Therefore, the simulations need to start by generating the exit products (protons), based on the energy and the spatial distribution of the primary reaction-inducing particles (neutrons).
For each separate excited state x in the daughter nucleus ( 12 B; see table 1) the neutron energy E is sampled from some preselected energy distributionφ E (x, E), where we use the hat-notation· to indicate the simulated (as opposed to the later determined, experimental) quantities. These distributions are best selected as uniform or isolethargic, for the simplicity of later analysis. The produced proton direction in the center-of-mass frame is then sampled from a preselected angular distributionÂ(x, E, χ), which is best selected as isotropic. The proton energy in the center-ofmass frame is calculated based on the Q-value for a particular excited state. The proton energy and direction are then boosted into the laboratory frame by the proper (in our case relativistic) transformations. As for the initial proton position, its radial distribution relative to the direction of the neutron beam must be sampled according to the known neutron beam profile; alternatively, the data need to be properly reweighed according to the same beam profile during the later construction of the E E matrix. The sampling (or the later data reweighing) of the initial proton position along the neutron beam direction depends on the properties of the simulated sample and may vary between extremely simple and rather involved. In case of the thin sample -implying the combination of the geometric thickness and the total cross section such that ηΣ tot (E) 1, as discussed in a context of eq. (2.6) -the longitudinal proton distribution may be sampled uniformly, as the neutron beam attenuation along the sample is negligible. This was the case with our setup. Otherwise, if the beam losses are known to be considerable, then the relative proton production probability along the sample must be properly accounted for. In case of the homogeneous and geometrically regular sample this correction amounts to the factor 1 − e −ηΣ tot (E)×d/D with d as the proton production depth and D as the sample thickness along the neutron beam direction; however, this procedure still does not take into account the multiple scattering effects. For more complex samples the correction involves its full spatial characterization.
Each coincidental proton detection by any pair of ∆E-E strips needs to be recorded by outputting the relevant physical parameters of that particular event. The necessary data consist of: (1) the primary neutron energy E; (2) the proton emission angle χ from the center-of-mass frame; (3) the unique designation of the activated ∆E-E pair of strips; (4) the energy deposited in those strips. In addition, the excited state x, the sampled neutron energy distributionφ E (x, E) and the proton angular distributionÂ(x, E, χ), together with the total numberN E (x) of generated protons within the sampled neutron energy interval E also have to be documented for a complete and meaningful utilization of the simulated data. By the virtue of eq. (2.2), the elements of the design matrix E E from eq. (2.28) may be treated as the integrals over the detected counts, so that by identifying the amount d 2N (x, E, χ) of generated protons as: where the branching ratios ρ(x, E) and angular distributions A(x, E, χ) are now taken to be known from an independent source of information.
We remind that the energy deposition cuts used in the analysis of the experimental data are to be implemented precisely at this point, in the construction of the matrix E E or the vector ì (E) , thus directly affecting the numbersN (E) αx of the acceptable counts. It is also worth noting that the weighting function w(E) is determined by the actual experimental conditions, as opposed to the arbitrary sampling distributionsφ E (x, E q ) andÂ(x, E q , χ q ). In that, it is evident that both the simulations and the computational procedures from Eqs. (A.3) and (A.5) are immensely simplified when the uniform neutron energy distributionsφ E (x, E) = 1/|E| and the isotropic proton angular distributionsÂ(x, E, χ) = 1/2 are used.

B Method extension
We shortly comment on the possibility of the further method generalization that may be applicable under specific conditions, namely the high statistics and at least a partial separation of the excited states in the deposited energy spectra. For the silicon telescope consisting of ∆E and E-layers, the entire two-dimensional ∆E-E spectra may be considered in the most general case, as we do here. For simplicity, figure 7 illustrates the basic idea on the schematic example of the one-dimensional, e.g. (E + ∆E)-spectrum. Evidently, if the excited states are sufficiently far apart in energy (as determined by the detector resolution), the spectra shapes may serve as an additional source of information to be exploited. In this case one defines the differential coincidental detection probability: