Bug in FirstGeneration

[fzeiser]

Issue by Caronthir Friday Aug 30, 2019 at 13:50 GMT Originally opened as https://github.com/oslocyclotronlab/ompy/issues/32

---

When the first generation method was translated from Fortran to Python, the paper describing the method was used as the main guideline. However, the paper contains a notational quirk where the normalization n_ij does not describe a matrix element, but a vector element n_i derived from n_j. In the Python version, this was interpreted as a matrix element and implemented as such, but this gives a wrong normalization.

Furthermore, there are questions whether this normalization is correct. It may be that a simple row-normalization is necessary. This was implemented in a test version, but led to strange results. It hints at a mistake in the iteration where the weights are being set, as the resulting matrix seems to be off by one.

Ann Cecilie is convinced the first generation method should work perfectly on toy examples, so this should be the goal.

[fzeiser]

Comment by Caronthir Friday Aug 30, 2019 at 13:58 GMT

---

There are restrains put on the weights that might be unphysical. Whether or not this is in fact true must be explored. The restraint of the weights being non-negative leads to convergence to a wrong FG-matrix if the all generation matrix is already close the FG-matrix, as Ann Cecilie and Magne reportedly experienced some years ago with MAMA.

[fzeiser]

Comment by fzeiser Tuesday Sep 03, 2019 at 08:53 GMT

---

Just to be sure: Where in Guttormsen1987 to you find this matrix element / vector n_ij? I've seen in in rewritten versions, like Larsen2011 Eq. (2), and Eq. (3).

• @cecgresonant: Can we post some results from eg. mama, where one sees that it works perfect, and what perfect means? So both an input script and plots/numerical results?

[fzeiser]

Comment by Caronthir Tuesday Sep 03, 2019 at 10:48 GMT

---

In Guttormsen1987 it is written as a vector n_i, while in Larsen there is the n_ij notation, which was interpreted as matrix element and implemented as such, both in my version and Jørgen's version.

[fzeiser]

Comment by cecgresonant Tuesday Sep 03, 2019 at 10:51 GMT

---

Note the text (Larsen et al., 2011, p.3 top right column): "The singles-particle cross section is proportional to the number of events populating levels in a specific bin, and thus to the number of decay cascades from this bin. We denote the number of counts measured for bins i and j in the singles spectrum Si and Sj , respectively. " So i and j both refer to excitation-energy bins, in this case neighbouring bins.

[fzeiser]

Comment by cecgresonant Tuesday Sep 03, 2019 at 10:54 GMT

---

In Guttormsen (1987): Page 1, right column, Eq.(2) and the text below: "The coefficients ni of eq. (2) are determined in such a way that the area of each spectrum fi multiplied by ni, corresponds to the same number of cascades. Thus, for singles particle spectra of constant cross-section we have ni = 1."

[fzeiser]

Comment by cecgresonant Tuesday Sep 03, 2019 at 10:58 GMT

---

I will dig up the old stuff from the systematics paper (Larsen et al. (2011))

[fzeiser]

Comment by cecgresonant Tuesday Sep 03, 2019 at 12:05 GMT

---

I found my "exact" case, it converges after 5 iterations, and the method works as it should (using Mama version 7.5.2). I also found my smoothed cases with a Gaussian (symmetric) smoothing and a realistic smoothing (Gaussian on Ex and NaI-resolution on Eg), and will go through them ASAP. I will make a short write-up as soon as I get time to do it (I am alone with the kids this week so time is unfortunately limited :) )

[fzeiser]

Comment by fzeiser Tuesday Sep 03, 2019 at 14:42 GMT

---

I'm still not convinced that it's coded up all too wrong. But here something on the normalization factors.

I get lost when thinking and rethinking this. At some point of time we said that the n_ij factor might lead to something like a normalization for each Ex. However, I don't think so any longer.

From the notation used, I agree that

for singles particle spectra of constant cross-section we have ni = 1.

If we take an easy test case, see #24, and say that we have Where the true fg matrix looks like:

1	1	1
1	1
1

and the all generations spectrum would look like:

3	2.5	1
2	1
1

Say that the population cross-section is 100 for each bin. Then we get an all-generations spectrum like this:

300	250	100
200	100
100

As stated above, the population cross-section is constant. Thus n_i or n_ij, however we call it, is 1. So when we multiply n_i * f_i , where f_i is the spectrum at Ex bin i, we will get just the same matrix out again. There is not normalization, or "flatness".

[fzeiser]

Comment by fzeiser Tuesday Sep 03, 2019 at 14:44 GMT

---

I suggest that we wait until we get the test-case matrix from AC and run the FG method on it. Then we see whether there is a difference.

[fzeiser]

Comment by cecgresonant Tuesday Sep 03, 2019 at 17:50 GMT

---

I think this is much easier to discuss on the whiteboard... The point is to normalise in such a way that we do not over- or under-subtract the spectra below due to variations in the population cross section in the excitation-energy bins. Also, the point is to obtain a primary spectrum that has average multiplicity 1. If the original spectrum has average multiplicity 3, we need to subtract enough of the underlying spectra so that it ends up (in the final primary spectrum) to have multiplicity 1. I attach slides that Magne and me used for the mini-workshop this summer at MSU, see pages 26-31, in particular pages 27-28 miniworkshop70Ni_NSCL_Aug2019.pdf . I suggest we all take a look and then find some time to discuss together.

[fzeiser]

Comment by fzeiser Wednesday Sep 04, 2019 at 13:44 GMT

---

By the way: The all-generation matrix above is wrong. I shall upload some code to generate a "correct" all-generations from a first generations later.

[fzeiser]

Comment by cecgresonant Thursday Sep 12, 2019 at 11:47 GMT

---

Yes you are right, the all-generations matrix is not correct. Maybe try with the one I sent by email some time ago?

[fzeiser]

Comment by fzeiser Thursday Sep 12, 2019 at 14:45 GMT

---

I'll check. But after fixinf the all generations-matrix to a correct matrix (using the ag_from_fg method), also the ompy implementation seems to work on something trivial. I'll check on your case when I find the time to.

See here: The fg is a fg of ones. I'll have to plots the color bars, too, but the result is at least quite (though not 100%) good.

[fzeiser]

Comment by fzeiser Friday Sep 13, 2019 at 10:00 GMT

---

Another issue here: it should be checked that the trial fg matrix is actually a first generation matrix, i.e. that the calculated multiplicity is 1. It requires some thought to create this. The above/below spectrum of ones does not fullfill this! [1 0 0] [1 1 0] [1 1 1]

and Eg_mid = Ex_mid = [0.5, 1.5 , 2.5]

[fzeiser]

Comment by cecgresonant Friday Sep 13, 2019 at 12:58 GMT

---

Yes indeed!!! I did a similar mistake the first time I made an artificial (very simple) matrix for the systematic-errors paper. With a wrong f.g matrix and with a wrong (inconsistent) all-generations matrix, it is not possible to get the correct answer. One has to think carefully through the branchings and the corresponding number of counts in the generated matrices.