Bug in FirstGeneration

[ErlendLima]

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.

[ErlendLima]

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]

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?

[ErlendLima]

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.

[cecgresonant]

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.

[cecgresonant]

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."

[cecgresonant]

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

[cecgresonant]

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]

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]

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.

[cecgresonant]

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]

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.

[cecgresonant]

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

[fzeiser]

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]

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]

[cecgresonant]

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.

[fzeiser]

Added examples in following notebook: https://nbviewer.jupyter.org/github/oslocyclotronlab/ompy_further_examples/blob/master/firstgeneration/fg_tests.ipynb

Get primaries from all generation spectra and compare to the tru primaries:

• Spectrum generated (by dicebox(?)) from AC
• Other simple mock spectra

@cecgresonant: can you send me the true primaries of that example? Then I will put them in the plot, too.

@magneg and @cecgresonant: What do you think about the results?

[cecgresonant]

Great! we can also test the cases where I put on some detector resolution later. I will also look into your mock spectra when I get time.

To your questions regarding my self-made, artificial case:

1. Generating the "data": I just constructed this myself using MaMa, no DICEBOX or anything like that. All levels (level E1, E2, E3, see figs) are populated 100 times each.
2. Primaries: Level E1: Ex = 3.5 MeV, only one gamma ray Eg = 3.5 MeV, so 100 counts in this peak Level E2: Ex = 6 MeV, primaries are Eg = 2.5 MeV, branching 33% = 33 counts, and Eg = 6 MeV, branching 67%, so 67 counts. Eg = 3.5 MeV is obviously secondary (from E1) Level E3: Ex = 8 MeV, primaries are Eg = 2 MeV, branching 30% so 30 counts; Eg = 4.5 MeV, branching 20% so 20 counts, and Eg = 8 MeV, branching 50% so 50 counts. See attached figures from my lectures in FYS4515.

After 5 iterations with the f.g. method in MaMa, I get that the higher-generations leftovers are of the order of 10^-6 counts (practically zero). I attach a copy of how I ran the f.g. method in MaMa. I will also provide the f.g. matrix later - I cannot attach the MaMa file here but I will send it by email. fgmethod_exactcase.txt

[fzeiser]

Well' I can quite easily generate the "real" fg bu myself then. Anyhow, you can see that the fg method with it's default parameters gives a rather close approximation to the actual spectrum.

[cecgresonant]

Yes, but I am a bit worried about the leftover of the 3.5-MeV transition for the E2 level and the 3.5-MeV and 6-MeV transition for the E3 level, see attached fig. This should have been very close to zero, but it is actually visible with the present range on the y axis. How many iterations did you run?

[fzeiser]

It's the statistical multiplicity estimation (with default values ;P). If I change to "total", I get basically nothing. I'll commit the wole notebook tomorrow. (for now: observe that the peaks you were concerned about vanished):

[fzeiser]

It seems to work fine now. I copied the example to the further examples repo. See the jupyter notebook: https://nbviewer.jupyter.org/github/oslocyclotronlab/ompy_further_examples/blob/master/firstgeneration/fg_tests.ipynb